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ABSTRACT 


A  digital  computer  program,  originally  developed  by  George  E.  Fink, 
capable  of  determining  the  in-and/or  out-of-plane  vibration  frequencies 
of  a  single  plane  piping  system,  using  the  basic  transfer  method,  is 
modified  and  augmented  by  the  writer.   The  program  is  modified  to  use  the 
writer's  method  to  reduce  the  amount  of  calculation  and  is  augmented  to 
use  Marguerre  and  Uhrig's  modifying  frequency  determinant  method  and 
Pestel  and  Mahrenholtz' s  remainder  method  for  higher  natural  frequencies. 

An  analysis  is  made  of  difficulties  encountered  with  the  original 
transfer  method  and  utilizing  the  modified  methods. 

The  program  accepts  branched  systems  and  non-stiff  intermediate 
supports.   A  discussion  and  an  explanation  of  how  to  use  the  program  are 
included. 
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CHAPTER  I 
INTRODUCTION 

1-1  General  remarks. 

Today,  in  increasing  measure,  natural  frequencies  of  a  system  capa- 
ble of  linear  vibration  are  being  determined  with  the  help  of  the  method 
of  transfer  matrices.   This  method  is  well  explained  in  the  recent  book 
by  E.  Pestel  and  Leckie  [3].* 

This  method  which  has  been  known  for  more  than  ten  years,  proceeds 
with  particular  simplicity  especially  with  the  use  of  digital  computer. 
However,  it  meets  with  important  numerical  difficulties. 

Numerical  difficulties  were  observed  by  George  E.  Fink  in  his  U.  S. 
Naval  Postgraduate  School  Master's  thesis  "Vibration  analysis  of  piping 
systems"  [4] .  He  used  the  straight- forward  method  of  transfer  matrices 
for  calculation  of  the  natural  frequencies  of  vibrating  piping  systems. 

The  purpose  of  this  paper  is  to  deal  with  numerical  difficulties 
and  with  the  remedies  which  have  been  developed  and  also  to  investigate 
reduction  in  calculation. 

1-2  Scope  of  work. 

Numerical  difficulties  are  encountered  at  higher  natural  frequencies 
and  for  systems  having  very  stiff  exterior  springs.   In  examining  these 
difficulties  the  paper  by  Marguerre  and  Uhrig  was  encountered  [2],  This 
tends  to  explain  difficulties  both  with  the  original  transfer  method  and 
with  the  Pestel  and  Mahrenholtz's  modified  methods  [1].   Although, 
Marguerre  and  Uhrig  suggest  several  different  methods  for  dealing  with 

♦Numbers  in  brackets  refer  to  the  bibliography,  page  37. 


the  difficulties,  those  which  appear  to  be  certain  to  eliminate  the 
difficulties  are  so  far  outside  the  spirit  of  the  transfer  method  that 
the  advantage  of  this  method  would  be  lost.   Therefore,  it  was  decided 
to  proceed  with  the  transfer  method   to  get  as  much  good  from  it  as  one 
could. 

The  writer  set  out  to  investigate  the  remedies  for  failure  at  high 
frequency.   The  digital  computer  program  "VIPIPE"  presented  in  the  Ap- 
pendices was  used  to  determine  the  natural  frequencies  of  in-and  out-of- 
plane  vibration  of  a  planar  piping  system  based  on  the  transfer  method. 
This  program  was  originally  developed  by  Fink  [4],  but  has  been  modified 
and  augmented  by  the  writer. 

1-3  Assumptions  and  limitations. 

1.  System  is  conservative. 

2.  Distributed  mass  or/and  lumped  mass  model  is  used.  For  the  lump- 
ed mass  model  used  in  curved  section  (or  in  straight  section)  sufficient 
subsectioning  is  carried  out  to  give  acceptable  answer. 

3.  The  transfer  method  surely  fails  if  the  system  has  rigid  exterior 
elements.   The  hangers  are  introduced  as  a  support  by  a  point  matrix  with 
appreciable  linear  and/or  torsional  compliance. 

4.  Program  VIPIPE  is  limited  to  a  system  of  a  maximum  of  fifty  com- 
ponents, twenty  branches,  and  twenty  hangers.   It  can  only  determine  in- 
and/or  out-of-plane  natural  frequencies  of  planar  system. 

5.  Boundary  conditions  should  be  such  that  half  of  the  state  quan- 
tities are  vanishing.* 

*Refer  to  Appendix  B. 


1-4  Notation. 

[    ]  matrix. 

■f   J  column  matrix. 

X,Y,Z  cartesian  right  handed  coordinate  system. 

u,v,w  displacement  in  the  XYZ  direction  respectively. 

g  -f  fy  rotation  about  XYZ  axis  respectively. 

M  bending  moment. 

N  normal  force. 

T  torque . 

V  shear  force. 

Z.  state  vector  at  location  i. 

x 

fZ.]  state  matrix  at  location  i. 
1 

[U.  J  transfer  matrix  from  i  to  i  +  1. 

A(U»  frequency  determinant. 

Z^n(l>0)  modified  determinant.* 

[N]  non-vanishing  state  matrix  of  Z  . 

Uj^.  element  of[U]. 

[S]  spring  matrix. 

[S]  purged  frequency  determinant  in  the  P-M  method.* 

[I]  unit  matrix. 

[0]  null  matrix, 

d  displacement  vector, 

[d]  displacement  matrix, 

p  internal  force  vector, 

[p]  internal  force  matrix. 

*ln  the  text,  the  following  abbreviations  are  used.  M-D  method  stands 
for  Modified  Delta  method  and  P-M  method  for  Pestel  and  Mahrenholtz's 
remainder  method.   These  will  be  explained  in  what  follows. 


(G) 

(Rl),(R2},(R3) 

Vo 

RM 

X 

P0      »    Pl     '    P2 

Xq  ,  X,    ,  x^ 
Y0  ,    Y(     ,X|     »A. 

D,     ,    D^ 

K 

X 

coordinate  transfer  matrix. 

square  sub -mat rices. 

column  vector  of  every  elements  value  of  1.  (£'•  p*  *  /. 

remainder. 

correction  factor  vector. 

assumed  state  quantities  at  the  starting  boundary. 

correction  factors  in  the  P-M  method. 

purging  factors  in  the  M-D  method. 

eigenvalue. 

eigenvector. 


CHAPTER  II 
TRANSFER  METHODS  AND  NUMERICAL  DIFFICULTIES 


2-1  The  state  vector  and  transfer  matrix. 
A.   State  vector. 

The  state  of  vibration  at  location  i,  specified  by  state  quan- 
tities, is  described  by  the  state  vector  Z  .   This  vector  is  a  column 
matrix  whose  elements  completely  describe  the  instaneous  displacements, 
both  rectilinear  and  angular,  of  that  point  from  its  quiescent  position, 
as  well  as  the  corresponding  rectilinear  and  angular  forces  in  the  member 
at  the  same  point  and  at  the  same  instant,  i.e.,  those  forces,  if  a  cut- 
ting plane  were  passed  through  the  point  at  the  instant  of  interest,  would 
have  to  be  applied  to  each  cut  face  to  prevent  relative  motion  between 
them. 

In  the  case  of  planar  piping  system,  the  state  vector  has  six  ele- 
ments, three  pertaining  to  displacement,  and  three  pertaining  to  force. 

The  state  vectors,  for  the  in-  and  out-of -plane  cases  respectively, 
as  used  in  VTPIPE,  are: 

u  displacement  in  X  direction, 
v  displacement  in  Y  direction. 
Xf        slope  in  X-Y  plane  with  respect  to  X  asix. 


'Lp 


M?   moment  about  Z  axis. 

Vu    shear  in  Y  direction. 

N    normal  force  in  X  direction. 


(2-1-1) 


Z0P= 


S  twist  about  X  axis. 

T  torque  about  X  axis. 

-w  (negative  of)  deflection  in  Z  direction. 

if 


(2-1-2) 


slope  in  X-Y  plane  about  X  axis. 
My.   moment  about  Y  axis. 
[V-,        shear  in  Z  direction. 
Thus,  the  state  vector  in  this  application  is  of  order  6x1,  while  the 
system  state  matrix*  is  6x3  and  the  transfer  matrix  of  a  section  is  6x6; 
however,  in  developing  theory  we  will  use  the  more  general  case  2rxl,  2rxr, 
and  2rx2r  respectively.   A  local  XYZ  right  handed  cartesian  coordinate 
system  and  the  elements  of  the  state  vector  are  arranged  as  in  Fig.  2-1-1. 


Central  axis 


A  <Jt  *J=  < 


-JJi 


w 


ijr 


M* 


/, 


^ 


^ 


Fig.  2-1-1 
B.   Transfer  matrix  and  transfer  method. 

From  the  state  vector  at  location  i  (node  i) ,  the  state 

*The  concept  of  "state  matrix"  is  introduced  in  Section  2-3. 
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vector  at  location  i  +  1  is  given  by 

zi.  +  l   "  Ktll'Zi.  (i  =  0,1,2,. ..,n)         (2-1-4) 

Here,  matrix  fUi+|]   relating  the  state  vectors  of  node  i  and  node  i  +  1 
is  termed  a  transfer  matrix. 

We  number  the  node  i  +  1  immediately  following  the  node  i  from  the 
left  end,  the  beginning  of  the  main  system.   From  the  recurrence  relation 
2-1-4,  one  easily  obtains 

Zn   '  (U„J  W-  •  •  N  •  N  •  M  •  Zo 

-  (0)  .  Zo  (2-1-5) 

which  is  a  linear  relationship  between  the  state  quantities  at  location 
0  and  n  of  the  system  (boundaries  of  the  main  system) .  We  call  this 
method  the  transfer  method. 

Since  of  the  2r  state  quantities  at  each  boundary  r  quantities  are 
vanishing,  only  r  are  established  by  the  boundary  condition  at  location 
0,  and  r  homogeneous  equations  are  established  by  the  boundary  condition 
at  location  n.   The  natural  frequencies  of  the  system  are  given  by  the 
zeros  of  the  determinant  of  the  corresponding  r  homogeneous  equations. 
Since  the  application  of  interest  is  to  the  determination  of  the  natural 
frequencies  of  piping  systems,  it  would  be  instructive  to  illustrate  the 
origin  and  nature  of  frequency  dependence  in  transfer  matrices.  For  thfi 
end,  let  us  construct  a  point  matrix  for  the  in-plane  case  connecting  Z ^ 
with  Z*  (z£  and  Z\     are  state  vectors  of  right  side  and  left  side, 
respectively,  of  a  point  mass  m  at  node  i) . 

To  begin  we  assume  that  all  of  the  elements  of  the  state  vectors 
vary  with  time,  each  having  the  factor  cos  oOlt  . 

11 


Noting  that  the  deflection,  slope,  and  moment*  are  continuous  across 


the  concentrated  mass  m  ,  so  that 


uR  =uu 

JL         A. 


vl-v).   .    »■=< 


From  the   free  body  diagram  of  Fig.    2-1-2A,   B,    the   following   relations 
are  established. 


u(t)   =  u.costOt 
v(t)   =  v.coslOt 

V^  •   cos  w  t-V^    cos  [*)  t 
N^  •   cos  tO  t-Hj[  •  cos  ^  t 


,±J  ■  -u^costO  t 
^#  -  -vu)acos0t 


•""-^ 


7n 


dim 

6P^ 


(2-1-6) 
(2-1-7) 

(2-1-8) 
(2-1-9) 


and  substituting  Eqn.  2-1-6  into  Eqn.  2-1-8  and  Eqn.  2-1-7  into  Eqn. 
2-1-9  we  have 


(2-1-10) 
(2-1-11) 
Now,  the  relation  of  the  state  vectors  of  the  right  and  left  side  of 


I      »  V  -L    +  mv  tO 


N  J   =  N^  +  mu  tf 


the  point  mass  m  in  matrix  form  is: 
-  K 


-  W  •  *t 


i.e. , 


u 

• 

V 

V 

*v 

v* 

N 

LR 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

m(^ 

0 

0 

1 

0 

0 

0 

0 

0 

1 

V 

IT 


N 


JLL 


[Up]  is  a  point  matrix  for  in-plane  vibration.   It  neglects  any  rotational 
inertia  that  the  point  mass  may  have;  however,  rotational  inertia  could 

*The  moment  is  continuous  only  if  there  are  no  point  elements  having  finite 
rotational  inertia. 
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be  accounted  for  by  including  anothernnon^etro>  term  in  the  matrix. 


vi 


l   Q 


v? 


inufa 


V 


Nj     Ni 


^O 


N?    M 


miOV 


Fig.  2-1-2A  Fig.  2-1-2B 

The  derivation  of  transfer  matrices  is  treated  in  great  detail  in  refer- 
ence [3]  and  transfer  matrices  for  various  components  of  planar  piping 
system  are  incorporated  in  VIPIPE. 


2-2  Failure  of  transfer  method. 

For  higher  natural  frequencies  and  for  systems  having  very  stiff 
exterior  springs  numerical  difficulties  are  encountered.   In  these  cases, 
the  numerical  value  of  the  frequency  determinant  is  given  as  the  differ- 
ence between  large  numbers.   So  practically  the  natural  frequencies  can 
no  longer  be  determined. 

From  a  mathematical  point  of  view,  in  case  the  transfer  matrix  [u] 

k 

has  a  single  dominant  real  eigenvalue,  the  columns  of  (Uj  become  more 
and  more  parallel  as  k  increases,  and  approach  the  first  eigenvector  x 
[6], [2 J  .   In  this  case  the  transfer  matrix  degenerates. 


£ 


(2-2-1) 


[U]  x,  =  Kt  x, 

This  could  be  the  case  for  the  lumped  mass  system  as  the  number  of  lumped 
masses  increases.   This  was  observed  in  the  straight  section  lumped  mass 
system,  but  for  the  curved  lumped  mass  system  this  was  not  observed. 
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As  the  frequency  increases  (elements  of  some  transfer  matrices  are 
functions  of  frequency),  the  same  phenomenon  was  observed  in  the  distri- 
buted mass  straight  section  transfer  matrix.   For  example,  for  the  sample 
problem  of  Appendix  E-2-1, 

At  frequency  of  20  rad/sec, 

2 
the  largest  eigenvalue:   ^,  ■  .11484022x10 

second  largest  eigenvalue:  X_#?U  ■  .99980414  +  jO. 01979160 

transfer  matrix  (U)  is: 


.99980413 

0 

0 

0 

0 

•  90535567V"5* 

0 

.25105587V1 

. 25986435V3 

.22549093V"2 

.14255229\° 

0 

0 

.30840073\"1 

.25105587V1 

.  26651055  V"4 

.22549093V'2 

0 

0 

.47565399\4 

.30070198x6 

.25105587  V1 

. 25986435V3 

0 

0 

.56219618V 

.47565399V4 

.30840073  V"1 

.  25105589V1 

0 

-.43265609V 

0 

0 

0 

0 

.99980413 

*V  means  10  ,  etc. 

In  this  case,  that  is,  at  this  frequency,  we  can  see  that  there  is  no 
phenomenon  of  the  columns  becoming  parallel.   As  the  frequency  increases 
this  single  dominant  real  eigenvalue  keeps  increasing,  and  at  frequency 
of  1700  rad/sec. 

the  largest  eigenvalue:  \=   .5937220  x  10 

second  largest  eigenvalue:  ,\2,  \      -  -.11136847  +  j. 99377918 

transfer  matrix  [u]  is: 
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-.11136847   1  0 

0 

0 

0 

.53482189X"5 

0         !  .14843074\10 

.13191200X11 

.12023306\5 

.10685242x6 

0 

0 

.  16701804\9 

. 14843074  \10 

.13528929\4 

.12023306\5 

0 

0 

.18324149\15 

.16284870\16 

.  14843074\10 

. 13 191 200\1 l 

0 

0 

.20618799\14 

. 18324 149\15 

.  16701804^ 

.14843074\10 

0 

- .  18465906\6 

o 

0 

0 

0 

.11136847 

Normalized  values  of  each  column  of  [u]  are: 


.6031032^7 

0 

0 

0 

0 

-.4802273x"5 

0 

.810028o\"6 

.8100279V6 

.8100280\"6 

.8100280\'6 

0 

0 

.91 1464 l\"7 

. 91146407 \"7 

.  911464  lV7 

.  9114641 V7 

0 

0 

1 

1 

1 

1 

0 

0 

.1125225V1 

.1125225V1 

.1125225\'1 

.  1125225V1 

0 

1 

0 

0 

0 

0 

1 

We  can  see  several  of  the  columns  are  almost  parallel  to  each  other.  From 
this,  it  is  easy  to  see  the  numerical  difficulties  we  should  encounter. 
This  depends  on  the  nature  of  the  transfer  matrix.  For  systems  made  of 
various  sections,  multiplications  of  different  transfer  matrices  may  result 
in  the  difficulties  taking  care  of  themselves. 

2-3  Reduction  in  calculation  in  matrix  multiplication. 

In  equation  2-1-5,  [u]  =   [Uj  >  [u^J  . . .  [U2]  .  [U,]   ,  so  we 
should  perform  matrix  multiplication  of  2rx2r  matrix  times  2rx2r  matrix 
successively. 

Now  using  the  fact  of  the  boundary  condition,  introduce  Z  as 
2rxr  matrix  called  the  "state  matrix.1'  This  is  formed  from  the  state 
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vector  as  follows,  all  elements  of  the  state  matrix  are  zero  except  that 
in  its  ith  column,  unity  appears  in  the  same  position  as  the  ith  general- 
ly nonzero  term  reading  down  the  state  vector.   Then  we  can  multiply 
2rx2r  matrices  by  2rxr  matrices  and  thus  can  reduce  the  calculation 
considerably. 

For  example,  in  the  in-plane  case,  for  a  built-in  left  end 


0 
0 
0 

1 

0 
0 


0 
0 
0 
0 
0 

1 


The  frequency  determinant  is  obtained  in  the  usual  way  from  the  boundary 
conditions  at  the  right  end. 

We  use  the  term  "Delta  method"  to  refer  to  this  method  which  was 
discovered  by  the  writer  during  the  course  of  this  study.   This  pro- 
cedure is  used  in  what  follows  instead  of  the  usual  procedure  of  the 
transfer  method. 


2-4  Pastel  and  Mahrenholtz's  remainder  method. 

Instead  of  the  frequency  determinant,  Pestel  and  Mahrenholtz  intro- 
duced a  remainder  [l] , [3]  .   This  method,  presenting  the  system  state 
matrix  as  2rxr  matrix,  reduces  the  amount  of  calculation  considerably 
compared  to  the  original  transfer  method,  and  claims  greater  accuracy  at 
higher  frequencies  because  the  remainder,  which  is  the  difference  between 
small  numbers,  is  itself  small.  This  will  be  explained  later. 
A.   Remainder. 

For  the  non-vanishing  state  quantities  of  the  state  vector  Z 
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which  depends  on  frequency  U)  we  use  estimated  values  which  approximate 


the  true  values,  following  which  the  Z   (i  ■  1,  2,  3,  ...,n)  are  deter- 
mined. 

In  the  state  vector  Z0 ,  since  for  eigenvibration  only  the  ratios  of 
state  quantities  are  of  interest,  one  of  the  non-vanishing  state  quanti- 
ties is  assigned  the  value  of  1. 

For  example,  for  the  in-plane  case,  at  certain  frequency,  a  canti- 
levered  pipe  which  is  free  at  the  left  end  and  fixed  at  the  right  end  has: 


o  o 

P,  +  X, 

o  o 

p,  +  x2 


(2-4-1') 


0 

0 

0 

(in  what  follows,  superscripts  0  and  1  denote  the  round  of  calculation; 
i.e.,  super-script  0  means  first  round  of  calculation  with  approximated 
state  quantities  in  ZQ  ,  and  superscript  1  means  second  round  of  calcula- 
tion with  improved  state  quantities  in  Zfl  ).   In  equation  2-4-1',  we 

assigned  the  value  of  1  to  the  non-vanishing  state  quantity  in  the  first 

o      o 
row  of  Z0,  and  the  quantities  of  P,   ,  P    are  the  previously  chosen 

approximations  and  X,   ,  X   are  correction  factors  which  are  unknowns. 

The  state  vector  Z    of  (2-4-1')  in  other  matrix  form  is: 


1 

0 

0 

■ 

1  ' 

*: 

1 

0 

< 

*; 

0 

1 

x° 

0 

0 

0 

0 

0 

0 

.    0 

0 

0   . 
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(2-4-1") 


For  the  actual  matrix  multiplication  Z   of  (2-4-1")  is  used.   After  a 


round  of  calculation,  the  approximated  state  quantities  P,   ,  Pa   will 
be  improved  in  the  form 


P°  + 


K  + 


xv 


(2-4-2) 


where  ?(,   and  %  are  corrections  obtained  by  this  process;  the 

difference  between  % ,   and  X(   ,  X^   and  Xa   will  be  shown  later. 
Then,  the  improved  state  vector  Z0  will  be: 

z:  =  fi 


(2-4-3) 

0 

0 

0 
When  we  start  first,  usually  P,   ,  P2   are  set  equal  to  1,  in  the 
absence  of  better  information.   From  the  same  method  as  (2-1-5)  using 


Zrt  of  (2-4-1"). 


[uj  .  [tvJ  ....[u2J  .  (u,J  .  z 


From  the  boundary  condition  at  location  it,  !.,#.,  from  the  three  vanishing 
elements  of  Z^  ,  we  get 

0  »  A(U>)'(N°)-  X°=  [SI'  X°         (2-4-4  ) 
[N]   :   nonvanishing  sub-matrix  of  [Z0]. 
X   :   correction  factor  vector. 
A(<\>)  :   frequency  determinant. 
[S]   :  A(u))-(N) 
For  the  cantilevered  pipe  of  the  example  in  the  in-plane  case,  the 
frequency  determinant  A(k))  ,  see  Ref.  [3]  ,  consists  of  the  top  left 
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3x3  submatrix  of  the  transfer  matrix  [U]  ■ 
Equation  2-4-4  may  be  written: 


0 


P°  0 


0 
0 


f  1 1 


1° 

;0 


^32 


U|3 

\x33) 


X 


U21       Ll22       Ua 
U»,        U&        U33J 

tin  tUiaP*  i-ul3p: 

I  tfei  +  ^p;  1 1^3  p; 

^2|    £«   5a3 
$31    632    $3> 
In  equation  2-4-5,  the  first  column  of  [s]  no  longer  contains  large 

numbers,  but  the  second  and  third  columns  remain  unchanged.  We  call 
this  "purging"  the  first  column.   From  equation  2-4-5  we  may  write, 


• 

*                > 

1 

x; 

/ 

\A 

(2-4-5) 


%\ 


$U         5 13 


Sja    ^J3 


8  "l 

! 
0 


(2-4-6) 


SS1       5»3? 

Equation  2-4-6  has  first,  three  equations  with  two  unknowns  (two  cor- 
rection factors) ,  and  second  it  has  the  purged  column  to  the  left  and 
unpurged  columns  as  coefficient  matrix. 
Equation  2-4-6  is  the  same  as; 


iV 

= 

'    Uq 

Ul3 

u» 

U-22 

u25 

U31 

u32 

U33 

p:+  x; 


(2-4-7) 


From  2-4-7  we  choose  any  two  equations.   If  we  choose  the  first  two, 
then  by  Cramer's  rule; 
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— 

Uii 

U.3 

— 

Uu 

Uii 

w° 

Uai 

lAz3 

•  -  p;  ■ 

-p; 

£- 

Uzi 

Ua| 

*i" 

U|2 

U, a 

Uiz 

U|3 

U:a 

Ujtf 

Uzz 

Ui3 

-U^-^M^, 


The   A  '8  are  defined  by  the  preceeding  expressions.   If  u3  =  Uj^, 
where  Oto  is  a  natural  frequency  of  the  system,  then  the  third  equa- 

o         o 

tion  will  be  satisfied  if  we  substitute  these  values  for  X  ,   and  X^  . 
Otherwise,  substituting  the  above  value  for  X,   ,  we  get,  from  the 
third  equation,  a  somewhat  different  value  for  X    ;  we  denote  this  by 
the  symbol  Y^ 

U31 


«_     U3aA2 


X 


U33A, 

We  define  the  remainder  R°(lO)   to  be  the  difference  between  these 
values. 


(2-4-9) 


r°(u>)=  Xi°-  y;  = 


A3         Uji  A2         U31 


U33A1 


U35 


(2-4-10) 


As  we  see  from  equation  2-4-10,  for  fixed  0  ,  the  remainder  is  indepen- 

-  o 

dent  of  the  assumed  state  quantities  P(   and  Pa   .   If  the  remainder  is 
zero,  u3  =  tOj^  since  all  the  prescribed  boundary  conditions  will  be 
satisfied  simultaneously. 

For  the  next  round  of  calculation  for  the  same  frequency,  form  the 
arithmetic  average  of  the  two  values  X  and  Y  for  the  correction  factor 

%    ,  i.e., 

0     0 

Then,  new  starting  values  are: 


K-  -j(K  +  r>) 
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p:-p;-*:- 


I      '  I  ■     'M      A\ 


o         o 

which  are  also  independent  of  P,   and  Pa  .    A  repeated  round  of  cal- 
culation gives 

x;=  x;=-p;--£  =  o 


therefore, 


2 


2 
Now,  we  see  that  the  iteration  is  completed  after  the  first  round  of 

calculation  and  the  remainder  is  no  longer  the  difference  between  very 

large  numbers,  because  the  remainder  is  of  the  same  order  of  magnitude 

i       i 
as  the  corrections  X^  and  Y2   .   During  the  iteration  seeking  the  natural 

frequencies,  we  calculate  R  (^\)     and   p^(  U)|)  's  (i»l,2),  and  increase 

frequency  by  an  amount  AO     to  get  lOi  =  0|+AO.   Then  get  R  (u)j)  using 

p;  l^a)  — P:  (^i)  (i  =  1»2).   The  above  example  is  for  the  case  r  =  3.   In 

other  cases  for  r  >  1,  the  spirit  is  the  same  (in  case  r  ■  1,  we  cannot 

use  the  P-M  method) . 

B.   Various  ways  to  get  remainder. 

There  are  many  ways  to  calculate  a  remainder.   Each  one  has  a 

different  value  at  the  same  frequency  (<£%ti/j£) ,  but  passes  through  zero 

at  0=w)£. 

Some  fail  and  some  give  good  results  (see  section  2-4-D) . 
The  possible  ways  arise  from  the  following  choices; 
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1.  Which  one  of  the  nonvanishing  state  quantities  will  be 
given  the  value  of  1  in  the  state  vector  ZQ. 

2.  Among  the  r  homogeneous  equations,  as  in  equation  2-4-5, 

o      o         o 
which  r-1  equations  will  be  chosen  to  get   X,   »  Xji   »•••*  r_|  . 

3.  Which  of  the  X°  's   (i  1,2,..., r-1)  in  the  rth  equation  is 
to  be  regarded  as  Y  and  solved  for  in  terms  of  the  other  X  's'. 

We  will  call  these  different  ways  by  the  name  "sub-procedures." 

C.   The  relation  between  a  remainder  and  the  frequency  determinant, 
From  the  equation  2-4-10,  the  remainder  is 

R(U)  )  .   a53A3-v  (132&2   -|-  Ub|A| 

From  equation  2-4-8,  we  see  that  A's  are  the  minors  of  an  element  of 
the  frequency  determinant 

where  M(u3j  )  means  the  minor  of  u3) 
The  remainder  is 

A(U)) 


R(0  )  = 


U33  A 


where  A(  0  )  is  seen  to  be  the  frequency  determinant. 

Expressing  a  remainder  in  general  in  terms  of  the  frequency  deter- 
minant ; 

R(U)  )  =  ,,A^  (2-4-12) 

Ufi   :  coefficient  element  corresponding  to  Y   of  rth  row. 

Afc,  :  minor  of  uhs  in  A  (0). 

u^g  :   element  in  rth  row  and  the  column  corresponding  to  state 


quantity  of  1  in  ZQ. 
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D.   Failure  of  the  P-M  method. 

In  case  the  denominator  of  the  remainder  is  zero,  the  P-M  method 
certainly  fails.*   (Refer  to  equation  2-4-10  and  2-4-12).   When  the  deter- 
minant in  the  denominator  becomes  nearly  equal  to  zero  or  changes  sign 
passing  through  zero,  the  remainder  is  difficult  to  evaluate:   at  certain 
frequencies  the  remainder  may  show  an  infinite  or  jumping  phenomenon. 

However,  this  does  not  necessarily  show  up  at  the  same  frequency 
with  different  sub-procedures,  because  this  phenomenon  depends  solely 
on  the  nature  of  the  transfer  matrices. 

Figure  2-4-1  and  Figure  2-4-2  show  that  while  one  sub-procedure 
shows  an  inifinite  phenomenon  at  certain  frequency,  another  sub-pro- 
cedure does  not  show  the  same  phenomenon  at  this  frequency. 

Pestel  and  Mahrenholtz  illustrated  their  method  with  the  case  of  a 
homogeneous  beam  (2r  =4).   It  leads  to  very  good  results  for  the  bend- 
ing frequencies.   The  reason  is,  because  r  =  2,  there  is  no  phenomenon 
of  indefiniteness  causing  failure  mentioned  above  and  due  to  purging 
effect.   The  P-M  method  does  show  improvement  also  for  the  case  of  r  great- 
er than  two. 

For  sufficiently  high  frequencies,  the  P-M  method  fails.   The  reason 
is  as  follows,  if  the  columns  have  already  become  parallel,  purging  does 
not  do  any  good,  and  if  round-off  error  reaches  critical  point,  where 
remainder  is  no  longer  reliable,  this  method  fails. 


*When  this  is  the  case,  program  VIPIPE  skips  to  the  next  sub-procedure, 
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R<tf» 


tee  862  884 


^ (rad/sec) 


X  SCALE  =*°£-*£ 
Y  SCALE-   =5-0 IE"  07 


Fig.    2-4-1 


RCvO  ) 


X  SCALE.  -  2,0  E  f02 
Y  SCALE   =5.0  £-07 


oteaie  6!5~ 

Fig.    2-4-2 


^(rad-     ^c) 
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Figure  2-4-1  shows  fourth  and  fifth  mode  frequencies  of  System  5*  at 
frequency  826.88281950  rad/sec,  and  1126.83799669  rad/sec  (P-M  method- 
subprocedure  A  with  double  precision  arithmetic).  Figure  2-4-2  shows 
infinite  phenomenon  at  frequency  939.00298809  rad/sec  while  Figure  2-4-1 
does  not  show  this  phenomenon  at  this  Frequency.   And  Figure  2-4-2  has 
print  out-put  with  mode  frequency  826  88281915,  also  with  939.00298809 
rad/sec  (P-M  method- subprocedure  B  with  double  precision  arithmetic). 
From  Figure  2-4-2,  we  see  that  the  latter  is  a  false  mode  frequency  and 
the  former  does  not  show  upon  the  graph  due  to  the  infinite  phenomenon 
of  the  latter,  but  we  know  that  there  is  a  mode  frequency  at  826.88281915 
rad/sec  (this  also  can  be  checked  with  the  frequency  and  remainder  value 
print  output) . 

2-5  Modified  Delta  method. 

Instead  of  purging  just  one  column  vector  among  r  columns  (as  in  the 
P-M  method)  of  frequency  determinant,  Marguerre  and  Uhrig  showed  in  their 
paper  [2]  a  possible  way  of  extending  the  purging  concept  for  problems 
greater  than  4th  order. 

Through  exactly  the  same  method  as  Delta,  we  get  a  frequency  deter- 
minant. At  high  frequency,  because  the  columns  approach  parallelism, 
the  value  of  the  frequency  determinant  tends  toward  zero.   If  the  columns 
are  exactly  parallel  to  each  other,  then  of  course  all  methods  fail.   If 
they  become  nearly  parallel,  then  with  a  finite  capacity  of  handling 
numbers,  the  results  become  meaningless. 

Using  Maguerre  and  Uhrig' s  purging  concept,  with  the  help  of  one 
column  of  the  frequency  determinant,  obtain  r-1  purging  factors  to  purge 

*See  Appendix  E-3. 
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another  r-1  columns.*  For  the  case  of  sixth  order  (r-3)  as  used  in 
program  VIPIPE,  the  frequency  determinant  at  a  certain  frequency  is: 

A(U>)  =  U„   u,a   Ub 
U2I   Uaa   u03 

Purging  factors  got  with  the  help  of  the  third  column  of  A(iO)  are 


-H 


an 


D 


U^L  +  _Ui 


3  v  Un 


t 


Lba 


IA23 


U33 


(2-5-1) 


Modifying  the  non-vanishing  sub-matrix  of  the  starting  state  matrix  [z0]  , 
we  have 


[N]  = 


1 

0 

0 

0 

1 

0 

Dl 

"D2 

1 

(2-5-2) 


Then,  for  the  same  frequncy,  repeat  a  round  of  calculation  to  get  the 
modified  frequency  determinant. 


A*(0)  = 


(2-5-3) 


uH  D,Ui3   Utt  D,u,3   U* 

U>l  D,U33     *A»  ftt^     U33 
From  the  properties  of  determinants,  A  (^  )  is  the  same  as  A,yn(^). 
But  the  first  and  second  columns  of  A™  (0 )  are  not  tne  same  as  those 
of  A(£>)  any  more.   If  the  columns  of  the  original  frequency  deter- 
minant were  exactly  parallel,  the  modified  columns  (first  and  second 
colum««)  of  the  modified  frequency  determinant  will  turn  out  to  be  all 

*We  use  the  term  "Modified  Delta  method"  to  refer  to  this  method. 
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zeros.   In  case  the  columns  were  nearly  parallel,  the  first  and  second 
columns  will  turn  out  to  be  very  small  numbers,  I.e.,  we  will  have  purged 
two  columns.   So  we  can  expect  a  more  accurate  value  of  the  frequency 
determinant  and  can  go  further  to  higher  natural  frequencies. 

If  D.  ■  D_  •  0,  it  degenerates  to  original  frequency  determinant, 
so  it  does  not  do  any  good.  If  D  =  D  ,  it  likely  means  that  first  and 
second  columns  are  the  same.   If  this  is  the  case,  this  method  just  fails. 

In  some  cases,  especially  systems  of  many  sections,  the  frequency 
determinant  does  not  show  parallelism.  Experience  shows  that  this  method 
gives  good  results  in  most  problems.   However,  this  method  also  fails 
at  sufficiently  high  frequencies. 

2-6  Intermediate  conditions  -  Branched  system. 

The  effect  of  a  joint  point  can  be  introduced  into  the  calculation 
by  means  of  a  single  point  matrix. 

This  section  is  mainly  for  the  derivation  of  the  branch  point  matrix 
of  the  Modified  Delta  method  and  the  Pestel-Mahrenholtz  method. 

The  derivation  of  branch  point  matrix  for  the  straight  forward  trans- 
fer method  (and  thus  also  for  the  Delta  method)  is  fully  explained  in 
Ref.  [3]  and  [4]  ,  but  a  brief  discussion  will  be  given  here  for  the 
purpose  of  refreshing  the  reader's  knowledge  and  moreover  explaining  the 
procedure  for  the  M-D  method  and  the  P-M  method. 

A.   Delta  method. 

Discontinuities  are  introduced  in  the  forces  whereas  the  main 
member  deflection  on  each  side  of  joint  are  the  same,  and  identical  to 
those  of  the  branch  at  the  same  joint. 
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JL    JR 


/^ 


0   -mcuTi  member    /  7,B 

branch 


Fig.  2-6-1 
By  a  suitable  coordinate  transformation,  it  is  possible  to  express  the 
force  system  of  the  branch  in  terms  of  the  main  system  coordinates. 
A  free  body  diagram  for  the  forces  at  the  joint  point  is  shown 
below.* 


7L 


H 


■3* 


Fig.  2-2-2A 
The  forces  at  the  joint  in  matrix  form  are. 


Fig.  2-2-2B 


'M' 

= 

'M' 

+ 

M 

V 

V 

V 

M. 

JR 

;N- 

JL 

.Mj 

(2-6-1) 


TIB 


The  displacement  is  continuous  at  the  joint,  so  Z  TR   can  be  obtained  by 
simply  adding  the  force  components  of  branch  state  vector  to  Z  ^  . 
The  column  matrix  of  nonzero  component  at  the  starting  boundary  will  be 
symbolized  by  VQ,     Here  V0   =  [l,  1,1} 


*Symbol  /\  indicates  that  the  vector  elements  are  expressed  in  terms  of 
the  branch  coordinate  system. 
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^ne 


d 
P 


TIB 


R, 

*2 


V, 


(2-6-2) 


From  equation  2-6-2 

"I        - 

The  matrix  product  TrJ'Tr.  ]    is  called  a  spring  matrix.  To  transform 
the  coordinate  system  of  the  branch  into  main  member  coordinates  at  the 
joint,  we  use 


here; 


d7.B  =  (frJ  '  ^ 

COS  0  -5IN0  0 

SIM  0  CO50  0 

0  0  I 


(2-6-4) 

I    0     0 
0   COS  0  -5IM0 
0   5IW  0       U)S>Cp\ 


where  <p     is  a  branch  joining  angle;  refer  to  Appendix  B-5. 
then; 

=  (5)  •  d^  (2-6-5) 

(  (sj  =  £  (*J  .   (R2MRi)  ~*[G]]>      ts)is  als°  called  a  sPrin8 

matrix.) 

Recalling  the  relation  at  the  joint 

d7R     sd7L     "  d7!B 

P^K     "^n.   +P71B     "PTL     +    Is)'    dTlB 
we  get  ZJR     in  matrix  form 
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3* 


TL 


I   o] -  fd  T 
s   ij  [p 

or        =   (Up]   .  Xyu 
Here  (Up)   denotes  the  branch  transfer  point  matrix. 


(2-6-6) 


B.   Modified  Delta  method. 

(In  this  section,  the  use  of  an  asterisk  *  will  denote  a 
purged  matrix  or  vector.   Also  note  that  the  quantities  d  and  p,  with  or 
without  *  and  ss  ,  represent  3x3  matrices  whereas  in  the  preceding  sub- 
section they  represented  3x1  vectors.) 

After  we  get  the  purging  factors  from^(u)),  the  new  starting  matrix 
at  location  0  is 
(<)  - 


r 


o 
i 

-D2 


At  the   joint 


(O 


d* 
P* 


fZ^u*J 


•    (N°J 


U-L 


For  a  branch,    let  the  non-vanishing  state  matrix  of  the  new  starting 
state  matrix     [  Z  0B]    be 

(nb]»      r    l       o       o 

0         10 
-DIB     -D26      1 

Here  D,6  ,  D^   are  branch  purging  factors,  the  calculation  of  which 
will  be  explained  later.   At  the  joint 


[£»*]    " 


d 

I  P  J 


•[■»}■        f  Ri 


TIB 


1R2J 


•      [   nb]= 


R  * 
R2* 
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[R, ], [R-J, [R*-]  &  fR*„]  are  defined  by  this  expression. 

i     z      1        i. 

Just  as  in  equation  2-6-5, 

-  [6il-[Ra]-tHJ-tNd"![Ri]"![&.] 
■   1&.]-[RJ-[R.  ]"[&,] 

-  [Si 

The  spring  matrix  turns  out  to  be  the  same  as  in  the  Delta  method. 
At  the  joint 

»?R  '  (P.J-IN1-V.+  [P*K 

[Sj  is  independent  of[NB]  ,  so  p*   is  independent  of  [NB]  . 


For  Z^*  we  have 


Z  * 


zi-  ^v  z- 


I    0 
S    0 

We  see  that  branch  transfer  matrix  [up I  is  the  same  as  in  equation  2-6-6, 
Purging  factors  can  be  obtained  from  [R- ]  orfR_)#  The  phenomenon  of 
columns  approaching  parallel  could  occur  in  [R_]  and  [R_] ;  VIPIPE  com- 
putes purging  factors  from  [R_]  ,  because  the  writer  was  inclined  to 
think  that  getting  purging  factors  from  [R  J  would  do  more  good  in  in- 
version of  [Rj  .   However,  [R  ]  could  have  been  used. 

Purging  factors  D|&  and  D2B  are  obtained  in  the  same  way  as  in 
equation  2-5-1. 

C.   Pestel-Mahrenholtz  method. 

At  the  joint  the  equilibrium  condition  must  be  satisfied, 
and  the  deflections  should  be  continuous. 
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-  d  '  = 


TIB 


R 

dug 


At  location  OB  (boundary  of  a  branch) . 


>0B 


•  X 


o 
06 


OB  J 


Here,  [  f\J  °    is  the  nonvanishing  state  matrix  of  [  Z^J and 
branch  correction  vector, 
let 


Kl=  f  P 


o 
OB 


Pi 


a8 


0 

I 

0 


0 

0 


Xqb- 


ofi 
o 


X 

Xi6 


Kb 


is   a 


Then,    the  state  vector  at  locatior   OB 

is 

0 

0 


Z    = 


£  OB    XO0 


P°       Vc 

r  IB    Ac 


r  0 

15    Aq6 


X 
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Pi  X 


0 

OB 


Xij 


for  the  built-in,    in-plane  case 


At  the   joint 

* 

£*.- 

'd 

'•NJ-x* 

ii- 

d 

•X" 

"p' 

JU 

.P* 

31? 

^TlB    ~ 

r  ^  i 
d 

•  [Nil  •  X^  = 

'R, 

■tNy-c= 

•  Pj 

.Raj 

(2-6-7) 


ft 


OB 
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where  fR  *)  and[R  *]  are 

For[N°]  and[X°)  ,  refer  to  section  2-4.  Also  [Rj  , [R  ]  are  expressed  in 
the  same  way  as  the  preceding  sub-section  B. 

(s*)  =  (e2)-(RD-(R,^[(r,) 

-  [<y-(Rj-KHN:Bp(R,|'-M 

-  [fr2]-[Ra]-(R.  ]"'(fri) 

-  (S) 

The  spring  matrix  also  turns  out  to  be  the  same  as  in  the  Delta  method. 
From  the  requirement  for  continuous  deflections  at  the  joint  we  get 
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We  see  that  the  branch  transfer  matrix  (Up]  is  also  the  same  as  in 
equation  2-6-6.  From  equation  2-6-8,  the  branch  correction  vector  is 
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After  we  get  the  main  member  correction  vector  (X  ) ,  we  can  get  branch 
state  correction  factors  in  terms  of  X  values.   As  explained  in  section 
2-4,  iteration  is  performed  in  the  same  way. 
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2-7  Computer. 

In  FORTRAN  63,  floating  point  quantities  have  an  exponent  and  frac- 
tional part.   In  single  precision  arithmetic,  the  fractional  part  has 
36  bits  and  the  precision  is  approximately  10  decimal  digits.   In  double 
precision  arithmetic,  the  fractional  part  has  85  bits  and  the  precision 
is  approximately  25  decimal  digits.   In  both  cases  the  range  of  numbers 
is  from  10"308  to  10308. 

In  the  FORTRAN  63  version  of  VIPIPE  (as  described  later),  either 
single  precision  (for  conservation  of  computer  time)  or  double  precision 
(for  greater  computational  accuracy)  can  be  selected.   In  the  FORTRAN  60 
version,  only  single  precision  is  available. 

The  natural  frequencies  of  a  planar  piping  system  are  found  using 
program  VIPIPE.   Program  VIPIPE  has  been  tested  with  FORTRAN  60  and  63 
using  digital  computer  CDC  1604.   The  number  of  natural  frequencies  that 
may  be  found  using  VIPIPE  and  the  accuracy,  especially  at  high  frequency, 
depends  upon  the  floating  point  significant  figure  capacity.  And  the 
significant  figure  capacity  required  for  any  system  is  a  function  of 
the  characteristics  of  the  system  itself. 
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CHAPTER  III 
DISCUSSION 

3-1  Accuracy  of  methods  and  assurance  of  the  solution. 

Confidence  in  accuracy  of  methods  was  established  by  comparing 
frequencies  of  several  systems  for  which  these  frequencies  could  be 
obtained  by  other  means,  i.e.,  by  means  of  either  a  closed  form  analyti- 
cal solution  or  by  recourse  to  the  Ref.  [5)  . 

The  accuracy  is  very  good  for  the  lower  natural  frequencies  in  the 
values  got  using  the  Delta,  method,  the  Modified  Delta  method,  and  the 
Pestel-Mahrenholtz's  method. 

In  most  cases,  up  to  fourth  mode  frequencies,  the  difference  be- 
tween the  system  natural  frequencies  obtained  by  analytical  means  and 
those  obtained  using  VIPIPE  are  less  than  0.004%.  As  frequency  in- 
creases we  got  more  accurate  natural  frequency  values  from  the  M-D 
method  and  the  P-M  method  comparing  to  the  original  transfer  method 
(Delta  method) . 

In  the  sample  problems  reported  in  Appendix  E,  we  have  found  that 
for  straight,  simple  configurations  (where  the  phenomenon  of  the  columns 
becoming  parallel  is  pronounced)  the  P-M  and  M-D  methods  give  one  or  two 
more  natural  frequencies  than  can  be  obtained  by  the  Delta  method  or 
the  original  transfer  method.  However,  in  more  complicated  configura- 
tions it  seems  to  be  possible  to  get  three  or  four  more  natural  fre- 
quencies.   In  the  simpler  cases,  there  are  comparison  values  available 
from  ,,exact"  theory  so  that  the  accuracy  of  the  results  may  be  establish- 
ed. However,  in  the  more  complicated  cases,  such  comparison  values  are 
not  available.  However,  we  have  reasonable  assurance  of  the  integrity 
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of  these  results  since  we  have  been  able  to  get  the  same  values  from 
the  M-D  method  and  from  several  different  P-M  subprocedures.   Thus, 
using  these  methods,  we  seem  to  be  able  to  double  or  triple  the  fre- 
quency range  in  which  valid  results  may  be  obtained. 

Details  of  the  systems  analyzed  and  the  results  obtained  may  be 
found  in  Appendix  E. 

3-2  Limitations  and  some  suggestions. 

Non-stiff  hangers  may  be  introduced  into  the  system  by  a  point 
matrix  with  appreciable  linear  and/or  torsional  spring  compliances. 
But,  for  stiff  hangers  the  transfer  method  fails,  so  we  cannot  use 
program  VIPIPE  for  natural  frequencies  of  a  system  which  has  stiff 
hangers.   In  particular,  the  suggestion  of  Fink  [4]  that  rigid  hangers 
can  be  successfully  handled  by  use  of  idealized  very  stiff  exterior 
springs  proves  to  be  untenable. 

Marguerre  and  Uhrig  suggested  possible  way  to  overcome  the  difficulty 
relating  to  stiff  hangers  by  introducing  unknowns  [2]  .   Introducing  un- 
knowns  means  method  abandoning  the  transfer  concept,  so  we  cannot  use 
VIPIPE  directly  to  solve  such  a  problem.  However,  the  writer  thinks 
that  there  is  a  possibility  of  modifying  VIPIPE  to  solve  such  problems 
by  getting  state  quantities  at  each  node  and  assuming  values  for  the 
unknowns.  For  details  refer  to  reference  [2]  . 
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APPENDIX  A 
DESCRIPTION  OF  PROGRAM 

A-l  General  remarks 

Program  VIPIPE  is  a  CDC  FORTRAN  63  language  digital  computer  program 
designed  for  use  in  the  vibration  analysis  of  planar  piping  systems, 
originally  developed  by  George  E.  Fink  [4]  .   The  program  VIPIPE  is 
modified  and  augmented  by  the  writer  to  use  Delta  method  (same  results 
as  original  transfer  method,  but  with  less  amount  of  calculation),  Modi- 
fied Delta  method,  and  Pestel  and  Mahrenholtz' s  remainder  method  with 
double  precision  or  single  precision  arithmetic. 

In-  and  out-of -plane  undamped  natural  frequencies  may  be  found 
through  the  use  of  program  VIPIPE.   How  this  program  may  be  used  is 
fully  explained  in  Appendix  B.   Considerable  flexibility  with  regard  to 
the  mathematical  model  employed  is  provided  in  the  program.   For  in- 
stance, one  may  use  combination  of  a  distributed  mass  treatment  for 
straight  section  with  a  lumped  mass  treatment  for  curved  sections  (this 
method  gives  the  greatest  accuracy  coupled  with  least  machine  time)  or 
a  lumped  mass  for  all  sections.   Also  one  may  consider  or  neglect  the 
effect  of  shear  deflection  and  rotational  inertia,  singly  or  together. 

Codes  for  the  boundary  conditions  and  various  methods  are  inter- 
preted within  the  program. 

Problem  solutions  proceed  by  iterative  processes  with  radian  fre- 
quency the  independent  variable.   The  sign  of  the  frequency  determinant 
or  remainder  is  used  to  control  the  search  for  and  convergence  to  a 
system  natural  frequency. 
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Output  and  transfer  to  the  next  method/or  transfer  to  the  next  pro- 
blem or  to  the  end  is  provided  in  the  program  when  it  reaches  signifi- 
cant figure  limitation  or  previously  set  iteration  limit,  and  also  when 
a  method  fails. 

All  subroutines  are  provided  in  binary  deck  in  single  precision 
and  double  precision  arithmetic  in  the  FORTRAN  63  version.   Modifications 
of  the  program  may  be  accomplished  by  changing  FORTRAN  cards 

only  in  the  main  program;  the  subroutines  are  in  binary  to  conserve 
space  and  time. 

A-2  Program  structure 

The  program  consists  of  a  six  section  main  body  and  twenty-seven 
subroutines.   In  the  FORTRAN  63  library,  GRAPH  subroutine  is  not  avail- 
able.  So  GRAPH  subroutine  is  provided  in  binary,  and  built  into  the 
program.   Also,  in  the  FORTRAN  63  library,  hyperbolic  functions  are  not 
available,  so  these  functions  are  computed  in  exponential  form. 

The  function  of  each  section  of  the  main  body  and  of  each  subroutine 
is  given  below. 

A.   Main  body  sections  and  their  functions. 

INPUT         Controls  read-in  of  data,  allocates  storage  loca- 
tions for  arrays  and  subscripted  variables,  de- 
cides the  first  method  to  perform. 
INVARIANTS  Computes  and  assigns  appropriate  subscripts  to  such 

system  component  characteristics  as  moment  of  in- 
ertia and  radius  of  gyration,  computes  frequency 
increment  for  the  first  two  modes,  and  sets  the  solu- 
tion acceptability  criterion. 
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CONTROL 
OUT-OF -PLANE 


ITERATION 


CONTROL       Constructs  transfer  matrix  of  the  system  for  the  in- 
IN- PLANE 

plane  case  by  calling  for  various  subroutines  in  a 

sequence  governed  by  the  system  geometry,  the  mathe- 
matical model  specified,  and  the  desired  assumptions. 
Performs  the  same  function  as  CONTROL  IN-PLANE  but 
for  the  out-of -plane  case. 

Evaluates  the  frequency  determinant  (modified  fre- 
quency determinant  in  the  M-D  method)  or  the  re- 
mainder and  utilizes  its  sign  in  controlling  the 
iterative  search  for  mode  frequencies.   Initiates 
output  when  the  specified  frequency  range  has  been 
fully  explored,  the  required  number  of  modes  found, 
the  specified  iteration  number  has  been  reached,  or 
the  significant  figure  limit  of  computer  reached. 
Controls  output  format,  provides  graph  output,  seeks 
the  method  to  perform. 

B.   Subroutines  and  their  functions. 

SUBSEC        Subdivides  components  to  be  treated  by  lumping  mass 
on  the  basis  of  L/D  ratio  and  starting  frequency  or 
required  mode  numbers.  (L/D  ratio  for  curved  section 
means  ratio  of  radius  of  curvature  to  diameter  and 
included  angle  of  arc,  as  applicable.) 
Constructs  the  transfer  matrix  of  those  components 
treated  by  lumping  mass  for  only  one  point  mass  and 
massless  subsection. 

Constructs  the  system  state  matrix  (6x3)  by  successive 
premultiplication  of  the  partial  system  state  matrix 


OUTPUT 


MATMUL 


FINMAT 
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FINBRA 


STATEM 


STATEB 


BRCOR 


DELMA 


by  the  transfer  matrix  of  each  component  as  soon  as 
each  of  the  latter  is  constructed.  For  the  lumped 
mass  system,  gives  successive  premultiplication  by 
the  transfer  matrix  constructed  by  MATMUL  as  many 
times  as  the  number  of  lumped  masses. 
Constructs  the  state  matrix  of  a  branch  system  by 
successive  pre-multiplication  of  the  partial  branch 
state  matrix  by  the  transfer  matrix  of  each  component 
of  a  branch  as  soon  as  each  of  the  latter  is  construct- 
ed.  For  the  lumped  mass  system,  gives  successive  pre- 
multiplication by  the  transfer  matrix  constructed  by 
MATMUL  as  many  times  as  the  number  of  lumped  masses. 
Constructs  starting  state  matrix  from  the  starting 
boundary  condition  of  main  member,  for  each  method. 
Constructs  starting  state  matrix  of  a  branch  system 
from  the  far  end  boundary  condition  of  a  branch  for 
each  method.  Computes  the  purging  factors  of  the 
branch  in  the  M-D  method. 

Constructs  branch  correction  matrix  for  the  assumed 
branch  starting  state  matrix  at  a  certain  frequency 
in  the  P-M  method.  Constructs  branch  purging  matrix 
in  the  M-D  method. 

Constructs  frequency  determinant  in  the  Delta  method. 
Evaluates  the  purging  factors  from  the  frequency 
determinant  of  the  Delta  method  and  constructs  the 
modified  frequency  determinant  after  another  round 
of  calculation  in  the  M-D  method.   For  the  P-M 
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method,  computes  the  correction  factors  for  each  sub- 
procedure.   For  each  sub-procedure,  if  it  fails  by 
the  fact  that  denominator  is  zero,  the  rows  are 
changed  and  a  second  test  is  made.   It  also  provides 
transferring  to  the  next  method  or  sub-procedure  if 
it  fails. 
DISTM,  DISTMO,  SFIELD,  SFIELO,  CFIELD,  CFIELO,  POINT,  POINO,  RIGID,  RIGIO, 
STIFCO,  HANGER,  HANGEO,  BRANCH,  BRANCO,  STAVEC ,  STAVEO,  INVERT. 
For  the  above  subroutines  refer  to  reference  [4]  . 

A- 3   Program  nomenclature. 

B1,B2,B3       Correction  factors  for  assumed  state  quantities  of 

the  far  end  of  a  branch. 
BC  Discriminant.   Controls  boundary  correction  in  the  P-M 

and  M-D  method. 
BL  Number  of  the  last  iteration. 

D1,D2  Purging  factors  in  the  M-D  method.   Correction  factors 

in  the  P-M  method. 
DEL  Array  of  numerical  values  of  frequency  determinants 

and  array  of  remainders. 
DV  Numerical  value  of  a  determinant  in  denominator  of 

XI  and  X2. 
EDL  Array  of  numerical  values  of  frequency  determinants 

or  remainders  in  single  precision  for  input  to  GRAPH. 
ERM  A  discriminant  same  as  REM  in  the  problem  statement 

data  card. 
FF  Discriminant.   Controls  output  to  provide  a  print 
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FINK 


HH 


HMH 

K 
KA 

KKK 


MAT 

NUMPTS 
NUMEND 
PI,  P2 
PMM 


output  that  analysis  is  terminanted  due  to  signifi- 
cant figure  limitation  in  the  P-M  method. 
Discriminant.   Controls  constructing  the  frequency 
determinant  or  modified  frequency  determinant  in  the 
M-D  method. 

Discriminant.  Controls  transferring  to  another  sub- 
procedure  or  to  another  way  of  testing  in  case  one 
fails  in  the  P-M  method. 

Discriminant.   Controls  subsectioning  in  lumped  mass 
system. 

Discriminant.   Controls  grid  in  the  GRAPH  output. 
Discriminant.   Terminates  iteration  when  specified 
iteration  number  is  finished. 

Discriminant.  Controls  transferring  to  the  output 
in  case  the  last  sub-procedure  in  the  P-M  method 
fails  or  controls  transferring  to  another  three  sub- 
procedures  (D,E,F)  after  three  sub-procedures  (A,B,C) 
are  finished  in  case  subprocedure  C  fails  in  Selec- 
tion 3  (for  Selection,  see  Appendix  B-9-6) . 
Discriminant.   Controls  matrix  multiplication  in  the 
lumped  mass  system. 

Number  of  points  in  one  graph  (except  the  last) . 
Number  of  points  in  the  last  graph. 

Assumed  starting  state  quantities  of  the  main  member. 
Discriminant.   Controls  three  methods  (Delta,  M-D, 
P-M  method)  specified. 
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PPM  Discriminant.   Assures  that  same  procedures  are  used 

for  out-of-plane  solution  that  were  used  for  in-plane 
solution. 

PM  Discriminant.   Controls  three  methods  to  be  used. 

PP  Discriminant.   Controls  the  three  sub -procedures 

(A,B,C)  in  the  P-M  method. 

RR  Discriminant.   Controls  testing  by  interchanging  two 

rows  in  case  one  sub-procedure  fails  because  the  de- 
nominator of  XI  or/and  X2  is  zero  in  the  P-M  method. 

R3  Array  name.   Matrix  product  [r  ]  .  [gJ  ,  used  in 

the  branch  correction. 

REM  Discriminant.  Controls  three  ways  of  selecting  sub- 

procedures  in  the  P-M  method. 

Q1,Q2,Q3       Assumed  starting  state  quantities  of  a  branch. 

SS  Discriminant.   Controls  the  three  sub-procedures  (D, 

E,F)  in  the  P-M  method. 

XX  Array  name.   Branch  correction  matrix  (3x3)  in  the 

P-M  method  and  branch  purging  matrix  in  the  M-D 
method. 

X1,X2,Y1,Y2    Correction  factors  for  assumed  starting  state  quanti- 
ties of  the  main  member  in  the  P-M  method. 

YY  Array  name.   Branch  correction  matrix  (3x3x23). 

! 
For  other  nomenclature  see  reference   [4] 
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APPENDIX  B 
INSTRUCTION  FOR  PROGRAM  USE 

B-l  General  remarks. 

Appendix  B  is  intended  to  provide  complete  instructions  for  using 
program  VIPIPE.  VIPIPE  is  FORTRAN  63  language  program  with  twenty-seven 
subroutines.  All  subroutines  are  provided  in  a  binary  deck  so  that  the 
compiling  time  is  reduced  as  is  the  volume  of  the  deck. 

Program  VIPIPE  can  be  used  in  single  precision  or  double  precision 
arithmetic.  Also  VIPIPE  is  provided  in  FORTRAN  60  language  with  single 
precision. 

The  Delta  method,  the  M-D  method,  and  the  P-M  method  are  all  incor- 
porated in  program  VIPIPE.   The  M-D  method  and  the  P-M  method  do  appear 
to  permit  obtaining  more  natural  frequencies  (one, two,  or  possibly  several 
more)  than  can  be  obtained  with  the  Delta  method  (same  result  as  the 
original  transfer  method).  However,  the  P-M  method  has  the  disadvantage 
of  introducing  an  odd  behavior  of  the  remainder.   In  some  cases,  this 
odd  behavior  can  easily  be  misinterpreted  so  as  to  lead  to  false  fre- 
quency determination.  However,  there  are  a  whole  group  of  associated 
P-M  modifications*  and  the  odd  behaviors  occur  at  different  frequencies 
for  the  different  modifications.  Thus,  by  making  the  calculation  more 
than  once,  using  a  different  modification  each  time,  the  true  natural 
frequencies  can  be  sorted  out. 

Program  VIPIPE  uses  the  Delta  method  as  the  basic  procedure.   If 
more  frequencies  are  required  than  have  been  obtained  by  this  method, 

*P-M  modifications  means  sub -procedures  of  the  P-M  method. 
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one  may  then  use  one  or  more  of  the  available  associated  P-M  modifica- 
tions and/or  the  M-D  method.   Graphical  output  will  assist  in  identifying 
the  natural  frequencies.   The  writer  recommends  to  use  the  P-M  modifica- 
tions as  original  method  in  getting  more  natural  frequencies,  and  to 
use  the  M-D  method  for  purposes  of  comparison. 

The  accuracy  and  number  of  natural  frequencies  that  can  be  obtained 
using  the  methods  depends  on  the  capacity  of  handling  numbers  accurately. 
As  we  can  see  in  Appendix  E,  double  precision  arithmetic  gives  better 
results  especially  for  complicated  systems  for  which  many  matrix  multi- 
plication are  performed.   But  using  double  precision  arithmetic  results 
in  considerable  increase  of  machine  time.   And  also  for  various  methods, 
the  machine  time  is  different.   The  Delta  method  and  the  P-M  method  each 
take  about  the  same  amount  of  machine  time,  but  the  M-D  method  takes 
almost  twice  this  time. 

So  a  user  should  consider  the  frequency  range  under  investigation  and 
the  nature  or  characteristics  of  the  system  as  well  as  the  machine  time. 

For  comparison  purposes  Section  B-9-6  of  this  Appendix  lists  machine 
times  for  one  iteration  for  some  example  problems  using  each  of  three 
methods  with  double  and  single  precision  arithmetic. 

Complicated  systems  are  divided  into  one  main  member  and  branches. 
System  components  are  defined  as  sections  of  pipe,  straight  or  of  con- 
stant curvature,  an  elbow,  a  spring  hanger,  a  coupling,  a  flange.   Dia- 
meters, wall  thickness,  density,  elastic  and  shear  modulus  may  vary  from 
component  to  component,  but  may  not  vary  within  a  component.   Both  tors- 
ional and  linear  spring  rates  (for  hangers)  must  have  finite  values. 
Thus  perfectly  rigid  hangers  cannot  be  treated. 
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B-2   Coordinate  system. 

A  local  coordinate  system  is  constructed  at  each  point  of  interest 
in  the  system.  A  right  handed  coordinate  system  is  used  whose  X  axis 
is  the  local  tangent  to  the  central  axis  of  the  main  member,  always 
directed  away  from  the  left  end.  The  X-Y  plane  coincides  with  the 
quiescent  plane  of  system.   The  Z  axis  is  always  directed  downward.  The 
direction  of  the  local  Y  axis  is  chosen  so  as  to  complete  the  right  hand- 
ed system. 


Fig.  B-2-1 
The  branch  coordinate  starts  from  the  end  of  the  branch  farthest  from 

the  joint  point  with  the  main  member.   The  plane  of  the  branch  coincides 

with  the  X-Y  plane.   The  positive  Z  axis  of  the  branch  and  main  member 

are  parallel  and  have  the  same  sense.  The  Y  axis  is  oriented  so  as  to 

complete  a  right  handed  system. 


B-3  Data  required. 

Extensive  properties. 

a.  outside  diameter  (D)  -  inches. 

b.  wall  thickness  (DI)  -  inches 

c.  radius  of  curvature  (RHO)  -  inches  (Negative,  if  center  of 
curvature  lies  on  negative  Y  axis). 
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d.  length  (TL)  -  inches  (for  straight  section  only;  for 
curved  sections,  length  is  computed  by  the  program). 

e.  arc  central  angle  (THETA)  -  degrees 

f .  linear  spring  constant  in  X  direction  (CLX)  -  pounds  per 
inch 

g.  linear  spring  constant  in  Y  direction  (CLY)  -  pounds  per 
inch 

h.    linear  spring  constant  in  Z  direction  (CLZ)  -  pounds  per 

inch 
i.    torsional  spring  constant  about  X  axis  (CTX)  -  in-lb/radian 
j.   torsional  spring  constant  about  Y  axis  (CTY)  -  in-lb/radian 
k.   torsional  spring  constant  about  Z  axis  (CTZ)  -  in-lb/radian 

Intensive  properties 

3 

a.  density  of  component  (AAMUj  -  lb/ft 

b.  elastic  modulus  (AE)  -  psi 

c.  shear  modulus  (AG)  -  psi 

B-4  Component  identification.* 

Components  are  identified  by  two  numbers. 

A.   All  components  have  a  component  sequence  number  starting  from 
the  most  left  component  of  the  main  member.   This  sequence  number  is 
identified  by  the  order  of  data  cards  containing  component  extensive  and 
intensive  properties.   This  can  be  understood  from  the  figure  B-4-1. 

*Aslo  see  pp.  B-5  to  B-7  of  reference   (4)  . 
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4.     -£ 


Fig.  B-4-1 
B.   A  branch  component  identification  number  (NNBR)  is  read  in 
as  data.   It  indicates  whether  it  is  in  a  branch  or  main  member,  if  it  is 
in  branch,  the  position  in  branch. 
NNBR  0    main  member 
NNBR  1    the  first  component  of  a  branch  (farthest  from  the 

intersection) 
NNBR  2    intermediate  component  of  a  branch 
NNBR  3    the  last  component  of  a  branch 
NNBR  -1   the  only  component  of  a  branch 
This  can  be  understood  from  the  figure  B-4-2 


Brwnch  2     ^ 

-i  / 


0 


.A 


I 


Fig.  B-4-2 
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B-5  Branch  joining  angle 

Branch  joining  angle  to  main  member  is  measured  counter-clockwise 
in  degrees  from  the  X  axis  of  the  main  member,  as  it  is  oriented  immedi- 
ately before  the  inter-section,  to  the  X  axis  of  the  branch  as  it  is 
oriented  at  the  point  of  intersection. 

B-6  Boundary  conditions* 

Five  common  boundary  conditions,  those  of  a  fixed,  free,  pinned, 
propped,  or  a  roller  supported  end  are  indicated  by  the  code  given. 

The  program  interprets  the  code  and  forms  the  required  state 
matrix  (6x3). 

Boundary  condition  code: 

a.  SV  1      fixed  end 

b.  SV  2      free  end 

c.  SV  3      pinned  end 

d.  SV  4      propped  end 

e.  SV  5      roller  supported  end 

f.  SV  0      other  than  one  of  the  above. 

(The  user  must  construct  the  necessary  state 
vector  and  cause  to  be  read  in  as  data) . 
This  code  is  used  for  both  in-plane  and  out-of -plane  cases. 

B-7  Control  cards. 

A  user  sets  up  for  execution  with  a  master  control,  a  FORTRAN  con- 
trol card,  and  combination  of  END,  FINIS,  EXECUTE,  transfer  cards,  bin- 
ary cards,  and  FORTRAN  cards  properly  placed  in  FORTRAN  63  deck  (JOB  and 
END  cards  are  used  in  FORTRAN  60).   Figure  B-7-1  shows  the  deck  assembly 

*For  further  details,  see  pp.  B-8  to  B-10  of  reference  [4). 
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for  FORTRAN  63  use. 


-&DATA 
-^DECK 


,  BINARY  DECK 
-t(SUBROUTlNE) 


Y FORTRAN  DECK 
(MAIN  PROGRAM) 


Fig.  B-7-1 


B-8  Data  deck  assembly. 

The  data  deck  is  ordered  as  follows: 

a.  problem  statement  card. 

b.  boundary  condition  data  card(s). 

c.  branch  joining  angle  data  card(s).   (deleted  when  system 
has  no  branches) 

d.  component  extensive  properties  data  card(s). 

e.  component  intensive  properties  data  card(s). 

The  user  may  obtain  solutions  for  any  number  of  problems  desired; 
in  this  case  each  problem  data  deck  is  assembled  just  behind  the  other. 
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Following  the  last  problem  deck  a  blank  card  should  be  added. 

B-9  Flexibility  of  program  VIPIPE. 

B-9-1  Double  precision  and  single  precision  arithmetic. 

Program  VIPIPE  can  use  double  precision  process  only  with  FORTRAN  63. 
Double  precision  is  used  in  system  state  matrix  (UU) ,  branch  system  state 
matrix  (W) ,  and  frequency  determinant  or  remainder. 

Type  double  is  declared  in  main  program  and  subroutines  FINMAT,  FIN- 
BRA,  STATEM,  STATEB,  BRANCH,  BRANCO,  BRCOR,  DELMA.   By  removing  cards 
declaring  type  double  the  program  can  easily  be  changed  into  single  pre- 
cision.  Subroutines  are  provided  in  binary  deck  with  single  precision 
and  double  precision  arithmetic,  so  with  corresponding  binary  subroutine 
deck  and  main  program  in  FORTRAN  with  or  without  a  card  (card  number 
000120)  declaring  type  double  we  can  use  single  or  double  precision  arith- 
metic. 

B-9-2  GRAPH. 

A  user  can  get  GRAPH  output  by  option.   The  values  of  the  frequency 
determinant  and  the  values  of  the  remainder  are  not  usually  of  the  same 
order  of  magnitude  throughout  the  frequency  range.  Moreover,  because 
of  the  infinite  phenomenon  of  the  P-M  method,  it  is  not  generally  useful 
or  convenient  to  get  one  graph  for  the  whole  frequency  range. 

1.  Number  of  points  in  one  graph  can  be  changed  by  changing 
NUMPTS  in  card  006090. 

2.  Can  get  graph  with  grid  or  without  grid  by  changing  card 

006100. 

K  -  1   with  grid 

K  ■  0   without  grid. 

52 


3.  Graphs  are  labeled  from  1,2,3,... up  to  10.  Graphs  beyond 
tenth  one  will  oe  labeled  as  10.  The  last  graph  is  label- 
ed LAST. 

4.  For  each  graph,  the  whole  frequency  range  is  moved  (except 
the  first  one)  so  that  the  first  point  of  each  graph  starts 
from  0.  The  corresponding  value  of  frequency  determinant 
(or  remainder)  for  the  first  point  and  NUMPTS  (NUMEND  for 
the  last  graph)  are  printed  in  the  print  output. 

B-9-3  Shear  area  form  factor. 

A  shear  area  form  factor  of  one  half  is  built  in  program.*  This  may 
be  changed  by  changing  card  000600. 

B-9-4  Starting  frequency,  frequency  increment  and  solution  acceptability 
criterion. 

Using  the  M-D  or  the  P-M  method,  we  are  expecting  to  get  more  natural 
frequencies  at  high  frequency.  Thus  choosing  the  proper  starting  fre- 
quency for  the  M-D  or  the  P-M  method,  following  the  Delta  method  can 
save  machine  time.   This  can  be  done  by  putting  the  wanted  starting  fre- 
quency in  problem  statement  data  card. 

To  compute  the  starting  frequency  increment,  the  program  constructs 
a  synthetic  straight  pipe,  equal  in  length  to  the  length  of  the  main 
member,  and  having  diameter,  wall  thickness,  and  intensive  properties 
equal  in  magnitude  to  a  weighted  average  of  those  properties  of  the 
composite  members.   One  quarter  of  the  fundamental  frequency  of  this 
synthetic  pipe  is  taken  as  the  starting  frequency  increment  (card  001390). 
And  one  thousandth  of  this  increment  is  taken  acceptability  criterion 
(card  001410).  When  the  iteration  reaches  this  solution  acceptability 

*See  reference  (7j  .  53 


criterion,  the  corresponding  natural  frequency  is  determined  using  linear 
interpolation. 

This  frequency  increment  is  for  the  first  two  mode  frequencies,  but 
the  acceptability  criterion  is  the  same  for  all  mode  frequencies.  After 
second  mode  frequency  is  found,  frequency  increment  is  one  tenth  of  the 
difference  between  the  last  two  mode  frequencies  found  (card  003950). 
These  can  be  altered  by  changing  cards.  The  writer  found  that  a 
smaller  frequency  increment  gives  better  results  especislly  at  high  fre- 
quency with  the  P-M  method. 

B-9-5  Subsectioning. 

Subsectioning  for  the  lumped  mass  treatment  of  components  is  carried 
on  in  subroutine  SUBSEC.  A  component  is  divided  into  maximum  of  twelve 
subsections. 

For  all  components  having  a  length  to  diameter  ratio  greater  than 
six,  the  number  of  subsections  is  computed  as  two  times  the  value  of  a 
parameter  HM  which  appears  in  SUBSEC  (HM  is  number  of  modes  to  be  sought) . 
If  HM  is  greater  than  six,  HM  is  given  the  value  of  six.   In  case  the 
starting  frequency  is  greater  than  800  radians  per  second,  HM  is  also 
given  the  value  of  six.  HM  can  be  controlled  by  a  parameter  HMH  in  the 
main  program  (card  001250).  If  HMH  is  zero,  HMH  does  not  affect  the 
value  of  HM.  By  changing  HMH  to  other  than  zero  we  can  change  the  value 
of  HM  equal  to  that  of  HMH. 

Poor  subsectioning  might  give  some  loss  of  accuracy  in  determining 
natural  frequencies,  especially  at  high  frequency.  Checking  the  natural 
frequencies  by  increasing  the  number  of  subsections  using  parameter  HMH 
is  recommended. 
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B-9-6  Machine  time  and  printing  line  limit. 

In  FORTRAN  63,  estimated  time  is  set  to  thirty  minutes  and  print- 
ing line  limit  is  set  at  30000  lines  in  VIPIPE.  These  limits  can  be 
altered  by  changing  the  master  control  card.   If  a  limit  is  reached, 
computation  stops  automatically. 

For  reference,  machine  time  for  one  iteration  of  some  example  problems 
is  given  below. 


Example  System* 

Delta  meth. 

M-D  meth. 

P-M  meth. 

Double 
Prec. 

system  1 

0.131  sec. 

0.247  sec. 

0.144  sec. 

system  4 

2.10  sec. 

4.22  sec. 

2.14  sec. 

system  5 

4.799  sec. 

6.90  sec. 

3.505  sec. 

Single 
Prec. 

system  1 

0.044  sec. 

0.117  sec. 

0.045  sec. 

system  5 

0.377  sec. 

0.968  sec. 

0.49  sec. 

B-9-7  Performing  various  methods. 

One  selects  the  methods  to  be  employed  by  use  of  the  problem  state- 
ment data  card.  Then  the  program  picks  the  first  method  in  the  INPUT 
section  of  the  main  program  with  a  discriminant  PM,  in  the  order:  Delta 
method,  M-D  method,  P-M  method. 

Finishing  a  method,  the  OUTPUT  section  of  the  main  program  opens  to 
the  next  method.  Finishing  all  specified  methods,  the  program  goes  to  the 
out-of -plane  portion,  to  the  next  problem,  or  to  the  end. 

In  the  P-M  method,  there  are  six  sub-procedures.   These  sub-pro- 
cedures are  controlled  by  discriminatns  PP,  SS. 

PP  :  0.   subprocedure  A,     SS  :  0.   sub -procedure  D, 
PP  :  1.   sub-procedure  B,    SS  :  1.   sub-procedure  E, 
PP  :  2.   sub-procedure  C,    SS  :  2.   sub-procedure  F. 
*For  example  systems  refer  to  Appendix  E. 
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Each  sub-procedure  has  another  possible  way  by  interchanging  rows  of  final 
determinant  (refer  to  section  2-4).   If  one  fails,  the  alternative  one 
is  tested.   If  one  succeeds,  then  go  to  the  next  sub -procedure.   There 
are  three  ways  of  selecting  sub -procedures  by  discriminant  REM. 

REM:  2.   three  sub -procedures  in  the  order:  A,B,C.  -  Selection  1 
REM:  1.   three  sub-procedures  in  the  order:  D,E,F.  -  Selection  2 
REM:  0.    six  sub-procedures  in  the  order:  A,B,C,D,E,F.  -  Selection  3 
If  one  wants  just  one  sub -procedure,  by  changing  PP  to  a  value  of  2 
(card  000710)  in  selection  1  only  sub-procedure  C  is  used,  or  changing 
SS  to  a  value  of  2  (card  001000)  in  selection  2  only  sub-procedure  F  is 
used.   If  one  wants  just  two  sub-procedures,  by  changing  PP  or  SS  to  a 
value  of  1  sub-procedure  B  and  C  or  sub-procedure  E  and  F,  respectively, 
are  used. 

B-9-8  Limiting  the  number  of  iterations. 

One  can  limit  the  number  of  iterations,  for  each  problem,  by  a  para- 
meter HL  in  the  problem  statement  data  card.   In  the  P-M  method,  because 
of  the  infinite  phenomenon,  the  frequency  increment  in  iteration  can  be- 
come very  small.   Thus  the  iteration  might  keep  going  on  indefinitely. 

The  program  permits  six  hundred  iterations  as  maximum.   The  writer 
found  three  hundred  iterations  to  be  adequate  to  get  nine  natural  fre- 
quencies with  the  acceptability  criterion  set  in  the  program. 

B-9-9  End  conditions. 

A  boundary  condition  state  vector  must  have  three  non-zero  and  three 
zero  elements.  Frequently,  however,  an  actual  condition  may  not  meet 
this  specification.   For  instance,  an  end  may  be  supported  or  connected 
to  a  piece  of  machinery.  Then  we  approximate  it  as  being  held  in  a 
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mounting  block  having  specified  linear  and  torsional  spring  rates,  and 
to  this  attach  a  fictious  massless  short  section. 

Through  the  use  of  this  sort  of  fictious  end,  a  proper  end  condi- 
tion state  vector  may  be  constructed  for  any  real  situation  with  little 
or  no  loss  of  accuracy. 

For  hangers  exhibiting  coupling,  systems  having  unusual  flexibility, 
and  flanges  and  couplings  refer  to  pp.  B-19  to  B-20  of  reference  (4)  . 

B-10  Format  of  data 
B-10-1  Problem  statement 

The  card  has  15  fields  under  control  of  the  following  field  specifi- 
cations. 

312,  4F3.0,  2F8.0,  4F2.0,  2F4.0 
Field  names  and  the  information  each  field  conveys  are  given  below  in 
order. 

A.   Field  names  and  content. 

number  of  components. 

number  of  branches. 

code  for  in-and/or  out-of -plane  solutions. 

number  of  modes  sought. 

code  calling  for  the  use  of  lumped  or  distri- 
buted mass  model 

6.  SD       code  calling  for  shear  deflection  inclusion  or 
neglect 

7.  RI       code  calling  for  rotational  inertia  effect  in- 
clusion or  neglect 

8.  OMEGAG    upper  limit  of  frequency  range  of  investigation 

9.  OMEGAI    lower  limit  of  frequency  range  of  investigation 
and  initial  iteration  frequency.   If  0  is  placed, 
it  will  begin  at  .01  radians  per  sec. 
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1. 

NS 

2. 

NBR 

3. 

IOP 

4. 

HM 

5. 

AMYK 

B. 


10.  GDABC 

11.  BRC 

12.  GRAPH 

13.  REM 

14.  PMM 

15.  HL 
Codes  used. 

1.  IOP  0 


2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 
10. 
11. 
12. 
13. 


IOP  1 
IOP  2 
AMYK  1. 
AMYK  0. 
SD  0. 
SD  1. 
RI  0. 
RI  1. 
GDABC  0. 
GDABC  1. 
GDABC  2. 
BRC  0. 


14.  BRC   1. 

15.  GRAPH  0. 

16.  GRAPH  1. 


code  calling  for  print  output. . 

code  calling  for  the  inclusion  or  exclusion  of 
branch  correction  factors  in  P-M  method  and 
purging  factors  in  M-D  method. 

code  specifying  whether  there  is  to  be  graph 
output. 

code  for  choosing  methods  in  P-M  method. 

code  calling  for  the  combination  of  3  methods  of 
Delta,  M-D,  P-M  method  which  a  user  wants. 

number  of  iterations. 


Both  in-and  out-of -plane  mode  frequency  sought 
(in  this  case,  in-plane  frequencies  are  sought 
first). 

In-plane  frequencies  are  sought. 

Out-of -plane  frequencies  sought. 

Lumped  mass  model. 

Distributed  mass  model. 

Consider  shear  deflection. 

Neglect  shear  deflection. 

Consider  rotational  inertia. 

Neglect  rotational  inertia. 

Print  solution  only. 

Print  solution  and  input  data. 

Print  solution  and  input  data  and  graph  data. 

Do  not  consider  branch  correction  factors  or 
purging  factors. 

Consider  the  branch  correction  factors  or  purg- 
ing factors. 

No  graph  output. 

Graph  output. 
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17.  REM  0.    All  six  sub-procedures  in  the  P-M  method  (A,B, 

C,D,E,F). 

18.  REM  1.  Three  sub -procedures  in  the  P-M  method  (D,E,F). 

19.  REM  2.  Three  sub -procedures  in  the  P-M  method  (A,B,C). 

20.  PMM  1.  Delta  method  only. 

21.  PMM  10.  M-D  method  only. 

22.  PMM  100.  P-M  method  only. 

23.  PMM  11.  Delta,  M-D  method. 

24.  PMM  101.  Delta,  P-M  method. 

25.  PMM  110.  M-D,  P-M  method. 

26.  PMM  111.  Delta,  M-D,  P-M  method. 

B-10-2   Boundary  conditions. 

Field  specification  22F7.3 

Boundary  conditions  are  read  in  separately  for  in-  and  out-of -plane 
cases.  When  both  cases  are  to  be  investigated,  in-plane  case  data  are 
read  in  first. 

N  +  2  fields  must  be  filled  in  (N  is  number  of  branches) .   The  first 
and  last  fields  are  first  and  last  boundary  condition  codes  of  the  main 
member.   Branch  boundary  conditions  are  entered  in  the  intermediate  fields 
in  accordance  with  the  branch  sequence  number.  When  a  boundary  condition 
is  not  describable  by  one  of  the  codes,  a  zero  must  be  entered  in  its 
data  field.   For  each  zero  an  additional  data  card  is  required.  The 
limitation  in  non-zero  state  quantitites  at  the  boundary  must  be  three 
(refer  to  B-9-9). 

B-10-3  Branch  joining  angle  (degrees). 
Field  specification  is:   20F7.3 
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B-10-4  Component  extensive  properties  and  branch  component  identifica- 
tion. 
A  data  card  for  each  system  component,  ordered  by  sequence  number. 
Field  specification  is:   3F5.2,  2F6.2,  6F7.2,  12 
Field  names  are:  in  order, 
D,  Bti  RHO,  TL,  THETA,  CLX,  CLY,  CLZ,  CTX,  CTY,  CTZ,  NNBR 

B-10-5  Intensive  properties. 

The  first  intensive  properties  data  card  has  four  fields  with  speci 
fication: 

F7.3,   2E10.2,   F2.0 
Field  names  are:  AAMU,  AE,  AG,  DD* 


*When  the  intensive  properties  of  all  components  are  the  same,  the  DD 
field  is  left  blank.  When  the  intensive  properties  are  different,  the 
DD  field  is  non-zero,  then  an  intensive  data  card  must  be  prepared 
for  each  of  the  remaining  components  (including  hangers).  Additional 
data  cards  are  identical  to  the  first  one  except  that  the  DD  field  is 
blank. 
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APPENDIX  C 
FLOW  DIAGRAMS 

C-l  General  remarks. 

Flow  diagrams  are  Included  for  all  parts  of  the  main  program,  and 
for  several  subroutines.   Since  CONTROL  OUT  OF  PLANE,  FINBRA  flow  dia- 
grams are  identical  in  structure  to  those  of  CONTROL  IN  PLANE  and  FINMAT 
respectively,  only  the  latter  are  included.   Reader  can  refer  to  the 
flow  diagram  of  STAVEC  and  STAVEO  of  reference  [4],   Since  the  flow 
diagram  structures  of  SUBSEC,  BRANCH  and  BRANCO  are  similar  to  those 
of  reference  [4],  they  are  not  included  here. 

Flow  diagrams  of  other  subroutines  that  are  not  included  in  this 
paper  are  relatively  simple  ones,  so  the  reader  can  refer  directly  to 
the  program. 

C-2  Index  of  flow  diagram. 
Main  program. 

Flow  Diagram.  Page. 

INPUT  63 

INVARIANTS  65 

CONTROL  IN  PLANE  66 

ITERATION  68 

OUTPUT  70 
Subroutines 

FINMAT  73 

STATEM  74 

STATEB  75 

BRCOR  76 

DELMA  77 
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C-3     Flow  Diagrams 
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INPUT   FLOW   DIAGRAM 


INM--OI 


IEAD:    Dimension 
statement 


READ:   Coded  case2 
^boundary  conditions 


READ:   Coded  state; 
vector  SV0P(l,j)  / 

I«1.6 


READ:    Coded  case 
boundary  condi* 


2£L 


READ:   Coded  state 
\vector  SVIP(l,j) 


READ:   Branch  joining 
angle  PHl(l)      I«1.NBR 


READ:    D(l),Dl(l),RHO(l),TL(l),THETA(l)] 

clx(i),cly(i),ctz(I),ctx(i),cty(i), 

CYZ(l),NNBR(l)        I«1,HS_ 


\READ:    AAMU,AE,/ 
\  AC  DP        / 


-J , 

READ:      AMU(l),E(l),G(l)| 
I~2,NS J 


0MEGAI=-.01 


INPUT   FLOW   DIACRAM(CONT.) 


I 


AMYK=AMYK-1. 

OMEGAO-OMEGAI 

FFAC«.5    . 

L«=l 

M  =  l 

AA-  BB  — CC  =  GG  =  EE—  0  . 

AAMU—  AL  -  AD  =-  AID  -  0 . 

KK=KKK=RR=  PP=KA=0 

NHN=NBR-M 


thetr(i)*tb:eta(i) 
dt(i)*-  Dl(l) 

BAMU(I)-AMU(I) 
DI(I)«DI(I)-2.0*DI(I) 
I«1,NS 


GO     TO   INVARIANTS 


ea- 


I  VARIANTS  FLOW  DIAGRAM 


FROM   INPUT 


X 


_& 


/clx(i)-cly(i)=clz(i)\ 
ctx(i)=«cty(i)=  ctz(i) 

I=-1,NS 


AMO(l)a*AMU(l)»l&*(D(l)»»2-Dl(l)»«  2 ) /( 6912.*32. 17*12.) 


THETA(l)=THETA(l)**/l80. 

tl(i)-theta(i)»absf(rho(i)) 


AJT(l)«JL*(D(l)**4Dl(l)**4)/32. 

AJ(l)=AJT(l)/2. 

R(I)-1./(AJ(I)*E<I» 

AIX(l)»SQRTP((l)(l)**2+Dl(l)**2)/8 
AIY(I)=»AIX(I)/1. 414214 


|z(i)*i._ 


o(ff^*0 
=0 


HHM=£M    PH=0. 
CALL   SUBSEC 


HM= 


HHM 


jinbr(i 
j*0 


=0 


AL»AL+TL(I) 

ad»ad+d(i)*tl(i) 
ee  =  ee+-e(i)*tl(i) 

^fAID=AID  +  Dl(l)*TL(l) 
AAMU—  AAMU  t~AMU(  I )  *TL(  I ) 


40 


j& 


i-i-t-i 


AD  -  AD/AL 

AID  =  AID/AL 

APR  =7^64.*(AD**4-AID**4) 

AD0M-FR/4.0 

DOM  -  ADOM 

CR»  FR*0. 00025 

PB«4.73**2*SQETP(EE*AFR/(AAM0*AL**4) 


BDOM  -  ADOM 
BR=CR 
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INVARIANTS   FLOV    PI AGRAM(CONT. ) 


? 


TIl(ll)=PHl(ll)*E/l80, 
II=1,NBR 


^OL 


MN  =11 


0  CONTROL  IN  PLANE 


21 


TO  CONTROL   OUT   OF  PLANE 


CONTROL      IN     PLANE      FLOW        DIAGRAM 


I    FROM   INVARIANT  Si 


P1^=P2=Q1=Q2=Q3=1 

D1  =  P2=-B1=B2  =B3=0 

FINK-  0. 


SVl(jfl)-SVIP(j,l) 
J»1.6 


=SL^T>^L 


SVE(  J,  1)  -SVIP(j,NMBR)|< — ^<|v(NMBR)^>^Q- 

J-1.6 


€Hp> 


i_o 


^cam  5553 


^CALL 


STAVEC 


-fCALL 


ST^VECJ 
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CONTROL  IN  PLANE  DIAGRAM  (CONT.) 


X  [  P(lWCLX(l)-»-~CLY(l)-t-CLZ(l)|^ 


*0 


CALL  RIQID1      |CALL  DISTHI  ICALL  STIFCOl 


I-I-H 
MAT=  0 
BC  =  Q, 


-»)CALL   HANGER 


ALL  STATEB 


CALL  STATEM 
1 


MAT«cO 
BC*1. 


ALL  CFIELDh 


I 


CALL  SFIELD 


ICALL  POINT  1 

rCALLMATMUU 
I 


^— »fCALL  FINBRAf X^BRO>> 


£  0 


<0 

I  CALL     IrANCHK- 


ITERATION    FLOW   PI APR AM 


6.0 


Numeric p1  v°lue  of 
frequency  determinant 


9 


GG-2. 

HDEL-DEL(L) 

OMEGA(lrtl)-  0MEG*.(l)  +  2.*CR 

"   A 


>® 


GO  TO:  Address 
6001  of  OUTPUT 


A 

KK-1 
T  M-M-1 

>/ 

GO   TO   OUTPUT 

DEL(L)=X2-Y2 
Dl«"  XI 
D2=(X2-f-Y2)/2. 


V- 


DEL(L)«X1-Y1 

Dl-(Xl"fYl)/2. 

D2-=-X2 


[TQ,guTPgT|  \0mMnl^QKEcI(i^) 


mQ 


*GG=-1. 


-•V- 


*)MEGA(  LT1 )  a  0MEGA(  L)  -C  R 

"^- 


omega(l)-omegao> Hpdel-del(l) 

y <-        =r()       ' 1 


>0 
£0 


PDEL-  DELCl 


>0       fPDEL     DEL("L-l)] 


^w-^^ggjiaiiiy1 


*<f- 


& 


ITERATION    FLOW   DIAGRAM(C0NT. ) 


DELC-DEL(L)  "1        ,       >Q 


BB=6.  ]f — ^-<PDE  L*  DE 


1  OMEGAC  sQMEGA(L) 

i <0  X  >Q_ 

(ajc-TohegaC  l-mT==  omegaClKdom  |<_J!1^omega(  lKomegag^hIm^m^ 


|OKEGA(ltl )  ~OMEGA(l)  -(  OMSGA(L) -0MEGAC)/2.k <BB> 


A 


(£)»— |CC=xlT|f-iiyOMEGA(L)-OMEGA(L-l)-CR) ^-^a) 


hnvnm(*\-  nmntr  ■    (OMEGA(l) -OMEGAC )»ABSF(PELC) 


OMEGAC -OMEGA(L-l) 
AA  =  1. 
BB=1. 
DELC=DEL(L-l) 


-/  omega(l)-omega(l-i)-crV— i 
sn x/ 


u^ 


BM—M 


m^-[m^  MtiV<fQ>— ^- 
M  i0 


>f 

Ito  out  put! 


>0 


DELC=0. 

AA=BB-CC=GG  =  0 
DOM  -  ADOM 


<U)ON=  (0MEGAM(M-1)-0MEGAM(M-2)  )*.T| 

[0MEGA(Lt-l)="0MEGAM(M-l)f-2.*CR 


GOi58oAddress 

CONTROL   IN  PLANE 


® — ^L»Ltl}-jBL 


OM EG  AO  ^OMEGA(liaj]-*(A) 


Ito  onTpnT|f \  ka  » 1 1 


®*— 


OMEGA(lrfl)  »  OMEGA(L)-.9*DOM 
DOMsr  .5* DOM 


GO  TO  Addressl002 
CONTROL  OUT  OF  PLANE 


6? 


OUTPUT  FLOW  DIAGRAM 


FROM   ITERATION 


I  PRINT;    Heading"]* — <™> 

i.. A  ,  i  \>0 

-ft  \S     >l 


PRINT:    IN-PLANE 
frequencies 


Ju 


PRINT:    OUT-OF-PLANE 
frequencies 

I 


PRINT:   Mathmatlcal 
model  used 


MH--M 
LM  =  L 


<0 


PRINT:  Delta 
method  used 


>0 


PRINT:    M-D 
method  used 


PRINT:    P-M  method, 
6   subprocedure  used 


PRINT:  Input  data 


* — 


>0 


'~L 


PRINT: Subsection  data 
graph  data, input  data 
= I 


Sort  OMEGA 
&  EDL 

GRAPH 


Address  6001 


omegam(m)~o. 

M=1,MH 


omega(l)=  0. 

L  =  ULM 
4. 


|delc=  omegac— 0. 


TO 


ss»ss-i.fg- 


2L 


ss-o. 

REM  =  -1. 


I 


<0 


OUTPUT   FLOW  -PI  AGRAM(C0NT..) 


d> 


=  0 


SS=»1. 


i 


MN— AA=BB=CC  =  GG  =  KK=  0. 
J   APOM—BBOM 
CR=BR 

OMEGA(l)=OMEGAI 
OMEGAO=OMEGAI 
DOM^iAPOM 
KKK=KA=FF=?=0. 
L=M=1 


® 


SS  =  2. 


©  © 


^PP~PP-1. 


-0     ^      >0 


PP  =  2. 


<o 


PP=1. 


<§ 


®<r 


SS«-1. 
PP  =  0. 

REM=-1. 


PP  =  0. 
REM=-1. 


I 


t> 


^0 


}  PRINT:  Failed  subprocedure 
in  the  P-M  method 


9 


OUTPUT   PLOW  _DIAGRAM(CONT.l 

IT 


->JREM-  ERM 


^MS>     i0 


SELECT:  Next 
specified  method 


A. 


•* 


/    Finished  sll       "A 
\specified  method      / 


YES 


°<mN  '° 


JJ-JJ-t-l 

1 

P(l)=Z(l>=  0. 
I      1,NS 

) 

f 

GO  TO:    Address  450  INPUT 

MN^O 
PP  =  0 
PM=PP« 


:    CONTROL  OUT  1 


'11 


SUBROUTINE  FINMVT  FLOW  DIAGRAM 


=  0 


>0 


V(J,K)=0. 
J»1.6  K»1.3 


R=Z-1. 


M  =  l 


V(j,KhV(j,KHu(j,L)*UU(L,K) 
J-1,6     K-1,3     1^1,6 ' 


V(J,K)  =  0. 
J   1,6  K  1,3 


V(j,KKv(j,K)tU(j,L)*UU(L,K) 
.    J=l,6     K~l,3     L-1,6 

I 


M-tt+1  f- 


UU(J,K)=:U(J,K) 
J-1,6  K-1.5 


<o 


>0 


V(J,K)»0. 
J»l,6  K-1,3 


*   _ 


V(J,K)  V(J,K)  a(j,l)*uu(l,k) 
J-1,6     K-1,3     W,6 


RETURN 


|  RETURN 
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SI^ROTJ^INE^^T.F^FLOW_DIA_GR_AM_ 


uu(j-,k)=o. 

J-1,5   K-1,3, 


<0 


P1-P1+D1 
P2»»P2-H)2 


|_w(j,jk)»i.J 


J  =  Jtl 


J  UU(Jt3)=P2h- 


uu(j;mHl. 


I — I  J»  Jfl 


UU(J,2)=D1 
'UU(J,3)=D2 


| RETURN  I 


74- 


SUBROUTINE  STATEB  FLOW  DIAGRAM 


Bl—  B2  —  B3  —  0. 
D1=D2=0. 


v(j,k)^yy(j,k,n) 

J=l,3  K  =  l,3 


<1— 


v(j,k)  =  yy(j,k,n) 

J-~l,3  K=l,3 


CALCULATE: 
Purging  factors 
Dl  &   D2 


W(J,2)«-D1 
W(j,3)»-D2 


* , 

CALCULATE: 

B1,B2,B3. 
For  subproc. 
A  &D 


CALCULATE: 
Bl,B2,B3. 
For  subproc, 
C  &F 


CALCULATE: 
B1.B2.B3. 
For  subproc. 
B  &   E 


I 


Q1~Q1*B1 
Q2=Q2*B1   B2f 
Q3S»Q3*B1    B3 

J- 


VV(J,K)»0. 


I 


& 
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SUBROUTINE   STATEB    ,  LOW  DIAGRAM(C0NT. ) 


$ 


Wfj,l)-Q3 

vv(j,m)-i. 


>0 


=  £> 


r(j,iH2 

r(J,M)=l. 


M  -  M  -t-1 


± 


>o   <^j=° 


RET 


40 


J?vN 


XO. 


til 


^M> 
>0 


W(j,l)=Ql]  [w(J,l)-l 


M=M1-  1 


*  J  =  J-t-1 


SUBROUTINE  BRCOR  FLOW  DIAGRAM 


=  0 


'-<$> 


*0 


J± 


xx(j,k)»w(j,k) 

J-1,3     K-1,3 


xx(i,k)=w(i,k) 

XX(2,K>»VV(3)K) 

XX(J,K>W(4,K) 

K-1.3 


RETURN 


>0 


JK 


UUU(l,K)=UU(l,K) 
UUU(2,K)=UU(3,K) 
UUU(3,K)=UU(4,K) 


=  6 


uuu(j,k>»uu(j,k) 

J-1,3  K-1,3 


XX(J,K)  =  0. 
J-1,3  K-1,3 


xx(j,k)*xx(j,kKR3(J,l)*uuu(l,k) 

J-1,3     K-1,3     L-1,3 


RETURN 


'tff" 


SUBROUTINE   DELMA  FLOW  DIAGRAM 


L_=X1_=  X2  =  Y1  =  Y2—DJ=*D2 -  0  | 


>0 


V(L,K)=UU(N,  K  )| 


CALCULATE: 
Purging 
factor  D1,D2 


>0 


KKK=1 


RETURN 


D(2,K)  =  V(3,K) 

D(3,K)  =  V(2,K) 

K^l,3 


RR-1. 
HH»1. 


Icalculate"~dv| 


RETURN 


=C 


|FINpi.l                |FINK-Q. 
I  RETURN  H 


D(J,K)«V(J,K) 
J»2.3  K~l,3 


* . 

CALCULATE  DV  \ 


=  £> 


j&O 


-it — 

CALCULATE 
XI  &   X2 


*0 


5 

77 


£ 

I  FP=a71 

IT 

[return] 


*r 


\^^0 


— >l 


lfo#* 


20 


it 


CALCALATE  YJ]       [CALCULATE  Yl~l 


rHBtuM 


RETURN 


.ss; 

-*0      V<0 


ss^iTITff^i. 


RR=0. 
HH=1. 


RETURN 


D(2,K)-  V(3,K) 

D(3,K)-.V(2fK) 

K=l,3 


CALCULATE  DV 


=  6 


£0_ 


<0 


D(J,K)  =  V(J,K) 
J-2,3  K=l,3 


CALCULATE  DV 


7 


RETURN 


CALCULATE 
XI  ft   X2 


^C 


-6 


t\  ^o 

1 

*6 

t 

' 

1 

|  CALCULATE 

Y2 

> 

t 

|  RETURN 

20 


RETURN 


F CALCULATE  Yl 


RETURN 


Id 


SUBROUTINE  DELM£   FLOW  DIAGRAMfCONT. ) 


=  t> 


*°<S 

b><* 

> 

t 

y^ 

f 

SS  =  2.      ^  =  2.1 

■ 

V 

RR  =  0. 
HH=1. 

► 

|  RETURN 

=0 


<D 


ytD 


*=D 


D(2,K)     V(3,K) 

D(3,K)      V(2,K) 

K     1,3 


D(J,K) 


'V(J.K) 


J^2,3     K=1.3 


|calculate"W1 


|CALCULATE~W 


RR=  1 . 
HH=1. 


i    [RE 


t1: 


RN 


FF-1. 


RETURN 


J£0_ 


£ 


*o 


CALCULATE 
XI  &  X2 


=0 


*& 


^0 


<o 


iO 


13 

CALCULATE  Y2J   [CALCULATE  Yl] 


RETURN 


RETURN 


7? 


APPENDIX  D 

PROGRAM  LISTING 

D-l  INDEX 

Main  Program. 

Section 

Sequence 

Number 

INPUT 

000000  - 

001070 

INVARIANTS 

001080  - 

001500 

CONTROL 

IN 

PLANE 

001510  - 

002450 

CONTROL 

OUT 

OF  PLANE     002460  - 

003390 

ITERATION 

003400  - 

004420 

OUTPUT 

004430  - 

007990 

Subroutines 

SUBSEC 

008000  - 

008520 

DISTM 

008530  - 

009090 

DISTMO 

009100  - 

009660 

RIGID 

009670  - 

009920 

RIGIO 

009930  - 

010180 

STIFCO 

010190  - 

010430 

STIFOO 

010440  - 

010680 

STAVEC 

010690  - 

010990 

STAVEO 

011000  - 

011260 

INVERT 

011270  - 

012000 

HANGER 

012010  - 

012180 

HANGEO 

012190  - 

012360 

CFIELO 

012370  - 

012880 

POINO 

012890  - 

013100 

Page 

82 

83 

84 

86 

88 

89 

96 

97 

98 

99 

99 
100 
100 
101 
101 
102 
103 
103 
104 
105 
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Section 

CFIELD 

POINT 

SFIELD 

SFIELO 

MATMUL 

FINBRA 

FINMAT 

BRANCH 

BRANCO 

STATEM 

STATEB 

BRCOR 

DELMA 


Sequence 
013110  - 
013570  - 
013800  - 
014160  - 
014430  - 
014570  - 
014860  - 
015150  - 
015790  - 
016480  - 
017030  - 
017790  - 
018110  - 


Number 
013560 
013790 
014150 
014420 
014560 
014850 
015140 
015780 
016470 
017020 
017780 
018100 
019500 


Page 

105 

106 

106 

107 

107 

108 

108 

109 

110 

111 

112 

113 

114 


D-2  Listing. 
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PROGRAM  VIPIPE 


INPUT 


DIMENSION  AJT(50) .AJ 
1DT(  50)  »D('50)  »DI  (50). 
2(50) »OMEGAM(20) »DEL< 
3LY(50).CLZ(50).CTX(5 
46.1).SV(22)  .SVIP<6.2 
520)  ♦P(20) .YYI3.3.23) 

TYPE  DOUBLE  VV.UU.V. 

NM  =  0 

JJ=1 
450  READ  1500.NS.N8R.IOP 
1REM.PMM.HL 

1500  FORMATI3I2.4F3.0.2F8 
IF  (NS)  600.600.601 

601  NMBR=NBR+2 

IF( IOP-1) 1501 .1501.1 

1501  READ  1502. ( SV ( I ).I=1 

1502  FORMATI22F3.0) 
DO  1505  J=1,NMBR 
IF(SV(J))1503. 1503.1 

1503  READ  1504,(SVIP( I.J) 

1504  FORMAT(6F2.0) 

1505  CONTINUE 

IF(IOP)  1506.1506.15 

1506  READ  1507, (SVO(  I),I  = 

1507  FORMAT(22F3.0) 
DO  1510  J-l.NMBR 
IF(SVO( J) )1508,1508. 

1508  READ  1509,(SVOP( I .J) 

1509  FORMATI6F2.0) 

1510  CONTINUE 

1523  IF(NBR)  1513.1513,15 

1511  READ  1512. ( PH I ( I ) .1= 

1512  FORMAT  (10F7.3) 
15130READ  1514. (D( I ).DI ( I 

1TX( I ).CTY( I ).CTZI I  ). 

1514  FORMAT(3F5.2.2F6.2»6 
READ  1515  .  AAMU.AE. 

1515  FORMAT  ( F7. 3  .2E10.2  . 
IF(DD)1516. 1516. 1518 

1516  DO  1517  1=1, NS 
AMU(.  I  )=AAMU 
E(I)=AE 

1517  G( I )=AG 

GO  TO  1520 

1518  AMU(1)=AAMU 
E(1)=AE 
G(1)«AG 

READ  1519. (AMU( I).E( 

1519  FORMAT  <F7. 3.2E10.2 ) 

1520  IF(OMEGAI)  1121.1121 


(50) .R(50) ,AIX(50) »AIY(50)»PHI (20I.SVI (6.1) » 
TL(50).Z<5O).G(50),E<50) »AMU( 50 ) .THETA ( 50 ) , RHO 
600) .OMEGA (600) »A(6.6).B(6.6).U(6.6).CLX(50),C 
0) ,CTY(50) ,CTZ(50),V(6,3)»VV(6,3),UU(6.3).SVE( 
2) »SVO<22) ,SVOP(6»22) ,NNBR(50) ,BAMU(50) ,THETR( 
,R3(3»3) ,XX(3,3),ITITLE( 12),EDL(600) 
X2,X1,Y2,Y1,P1,P2,D1,D2,DEL 


.HM.AMYK.SD.RI .OMEGAG .OMEGA  I ,GDABC ,BRC .GRAPH, 
.0.4F2.0.2F4.0) 


506 

,NMBR) 


505 
.1-1.6) 


23 
l.NMBR) 


1510 
•1-1.6) 


11 
1.N8R) 

) ,RHO( I ) »TL( I ) .THETA ( I ) .ClX( I ),CLY( I ) .CLZ( I ).C 

NNBR( I  )  .1-1. NS  ) 

F7.2.I2) 

AG.DD 

F2.0) 


I  )  .G(  I  1.U2.NS) 
.1521 


000000 
000010 
000020 
000030 
000040 
000050 
000060 
000070 
000080 
000090 
000100 
000110 
000120 
000130 
000140 
000150 
000160 
000170 
000180 
000190 
000200 
000210 
000220 
000230 
000240 
000250 
000260 
000270 
000280 
000290 
000300 
000310 
000320 
000330 
000340 
000350 
000360 
000370 
000380 
000390 
000400 
000410 
000420 
000430 
000440 
000450 
000460 
000470 
000480 
000490 
000500 
000510 
000520 
000530 
000540 
000550 
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1121 
1521 


6091 

6092 

6093 
6094 

6095 

6096 
6097 

6098 

6099 
7000 

7001 
7071 

7072 
7073 


1522 


\ 


OMEGAI=.01 

OMEGA ( 1)=0MFGAI 

AMYK=AMYK-1.0 

OMEGAO=OMEGAI 

FFAO.5 

L«=l 

MM 

MN=0 

AA=0.0 

BB=0.0 

CC=0.0 

FF*0.0 

GG=0.0 

RR*0.0 

KK*0 

PP-0. 

KA  =  0 

(CK(C«0 

AAMU-0.0 

EE-0.0 

AL*0.0 

AD=0.0 

NNN=NBR+1 

IF(PMM*0.1-1.)    6091*6092*6093 

PM=-1. 

GO    TO    7000 

PM=0. 

GO    TO    7000 

1F(PMM*0.01-1.)    6094.6095*6096 

PM=-1. 

GO    TO    7000 

PM«1. 

GO    TO    7000 

IF(PMM»0. 01-1.1)    6097.6098.6099 

PM=-1. 

GO    TO    7000 

PM=0. 

GO    TO    7000 

PM»-1. 

PPM=PM 

IF(PM)     7073.7073.7001 

1F(REM-1.)     7071.7072.7071 

SS*-1. 

GO    TO    7073 

SS  =  0. 

DO    1522    I-l.NS 

DT(I)*DI(I) 

BAMUU  )"AMU(I) 

THETR( I)»THETA( I) 

DH  I  )=0(1  >-2.0*DI(  I) 


INVARIANTS 


DO  108  IM.NS 

IF(CLX( I)+CLY( I )+CLZ( I )+CTX( I )+CTY< I )+CTZ( I) ) 108.101.108 


000  560 

000570 

000580 

000590 

000600 

000610 

000620 

000630 

000640 

000650 

000660 

000670 

000680 

000690 

000700 

000710 

000720 

000730 

000740 

000750 

000760 

000770 

000780 

000790 

000800 

000810 

000820 

000830 

000840 

000850 

000860 

000870 

000880 

000890 

000900 

000910 

000920 

000930 

000940 

000950 

000960 

000970 

000980 

000990 

001000 

001010 

001020 

001030 

001040 

001050 

001060 

001070 

001080 

001090 

001100 

001110 


S3 


101  AMU( I)=AMU(  I )»3. 141 592654* <D( I )»*2-DI ( I ) **2 ) / ( 69 12 . *32. 1 7* 12  .  ) 
IF(RHO( I ) ) 102. 103.102 

102  THETAI I )=THETA( I ) »3. 141592654/ 1 80. 
TL< I >=THETA( I )»ABSF(RH0( I ) ) 

10  3  AJT( I)=3.141592654»(D< I )**4-Dl< I >*«4)/32. 
AJ( I )=AJT( I )/2. 
R( I )  =  1.0/(AJ( I >#E<  I  )  ) 
A I X { I )=SQRTF( (D( I )»*2  +  DI ( I ) **2 ) /8.  > 
AIY( I )=AIX<  I  1/1.414214 
IF(AMYK)104,105.104 
Z(  I  )=1.0 

IF<RH0( I) ) 105.106.105 
HHM=HM 
HMH=0. 

CALL  SUBSEC  ( RHO( I ) . THETA < I ) .0(  I > . TL ( I ) .HM.Z (  I ) .OMEGA  I ,HMH) 
HM*HHM 

IF(NNBR( I ) ) 108,107.108 
AL=AL+TL( I) 
AD  =  AD+D( I >*TL<  I  ) 
EE  =  EE  +  E(  I  )*TL(  I  ) 
AID  =  AID+DI ( I )«TL( I  ) 
AAMU  =  AAMU+AMU( I >*TL(  I  ) 
CONTINUE 
AD=AD/AL 
AID=AI0/AL 

AFR=3.141592654/64.*(AD**4-AID**4) 
FR  =  4. 730040 19**2*SORTF(EE*AFR/(AAMU*AL**4)  ) 
ADOM=FR/4.0 
DOM-ADOM 
CR=ADOM*0.001 
BDOM=ADOM 
BR  =  CR 
ERM=REM 

DO  109  11=1, NBR 

PHI ( II )=PHI ( I  I )*3. 14 1592654/ 180.0 
IF(IOP-l)  110.999,998 
MN=1 


104 


105 


106 
107 


108 


109 


110 


CONTROL   IN  PLANE 


999  Pl=1.0 

P2=1.0 

D1  =  0. 

D2  =  0. 

01*1.0 

02=1.0 

03=1.0 

B1=0. 

B2  =  0. 

B3  =  0. 

FINK=0. 

IF(SV( 1))1011, 1021, 1011 
1011  CALL  STAVEC  ( SV ( 1 ) , 1 ,SVI ) 

GO  TO  1041 
1021  DO  1031  J=l,6 


001120 
001130 
001140 
001150 
001160 
001170 
001180 
001190 
001200 
001210 
001220 
001230 
001240 
001250 
001260 
001270 
001280 
001290 
001300 
001310 
001320 
001330 
001340 
001350 
001360 
001370 
001380 
001390 
001400 
001410 
001420 
001430 
001440 
001450 
001460 
001470 
001480 
001490 
001500 
001510 
001520 
001530 
001540 
001550 
001560 
001570 
001580 
001590 
001600 
001610 
001620 
001630 
001640 
001650 
001660 
001670 


H 


1031  SVI (J*1)*SVIP(J*1)  001680 

1041  IFISV(NMBR) )1051»1061»1051  001690 

1051  CALL  STAVECISVINMBR) »1.SVE)  001700 

GO  TO  1081  001710 

1061  DO  1071  J=l»6  001720 

1071  SVEI Jtl>=SVIP<J.NMBR)  001730 

1081  JK=0  001740 

IF(NBR)1000, 1000, 1085  001750 

1085  00  1087  N=2.NNN  001760 
IF(SV(N) )1087, 1087,1086  001770 

1086  CALL  STAVEC(SV(N),N,SVIP>  001780 

1087  CONTINUE  001790 
1000  N=2  001800 

11*1  001810 

HH«0.  001820 

DO  1001  1  =  1, NS  001830 

MAT=0  001840 

BC*0.  001850 

P(  I  )=CLX( I)+CLY(I )+CTZ(I )  001860 

IF(P(I ) )37, 38,37  001870 

37  CALL  HANGER(CLX( I ),CLY< I) ,CT2( I ),U)  001880 
GO  TO  62  001890 

38  IF(AMYK)39,41,39  001900 

39  IF(RHO< I) )40.42,40  001910 

40  IF(Z( I ) )43. 45,43  001920 

41  IF(RHO( I) )40, 48,40  001930 

42  IF(ZU)  )46,47,46  001940 

43  CALL  CFIELD  ( R( I ) ,Z ( I )  ,G( I ) ,SD,RHO( I ) ,THETA( I ) ,D( I ) ,DI ( II , E( I ) ,A)    001950 

44  CALL  POINT  ( AMLU  I ) ,A I Y ( I ) ,OMEGA( L ) ,TL( I ) ,Z( I ) ,R I ,B )  001960 
CALL  MATMUL  (A,B,U)  001970 
MAT=1  001980 
GO  TO  62  001990 

45  CALL  STIFCO  < THETA( I ) ,RHO( I ) »LM  002000 
GO  TO  62  002010 

460CALL  DISTM  ( AMU ( I ) , A  I Y ( I ) , OMEGA ( L ) ,TL ( I ) ,R< I > ,0  I -I ) ,DI ( 1 ) ,G< I ) ,S0,R   002020 

II, E(  I  )  .FFACU)  002030 

GO  TO  62  002040 

47  CALL  RIGID  (  AMUU  )  ,TL  (  I  )  , OMEGA  <  L  )  ,AI  Y  (  I  )  ,R  I  ,U  )  002050 
GO  TO  62  002060 

48  IF(Z(I ) )49,47,49  002070 

49  CALL  SFIELD  (  DI  (  I  )  »TL (  I  )  »R < I  1 ,Z ( I ) ,G< I ) , SD,D( I > ,E ( I ) ,FFAC, A)  002080 
GO  TO  44  002090 

62  IF(NNRRU))  81,80,81  002100 
80  IF( 1-1)63,79,63  002110 
79  CALL  STATFM(SVI,UU,D1,D2,P1,P2,PP»PM,SS)  002120 

GO  TO  63  002130 

63  CALL  FINMAT(U,UU,V,MAT,A,Z( I ) )  002140 
DO  64  J  =  l,6  002150 
DO  64  K»l,3  002160 

64  UU( J»K)=V< J,K)  002170 
IF(PM)  88,1113,84  002180 

1113  IF(FINK)  88,84,88  002190 

84  IF(BC)  88,88,1111  002200 

1111  IF(BRC)  88,88,89  002210 

89  CALL  BRCOR(R3»UU,XX,JK,VV»PM)  002220 

MMM=N-1  002230 


B5 


DO  90  J  =  l,3  0022*0 

DO  90  K=1.3  002250 

90  YY< J.K.MMM>=XX< J.K)  002260 

88  GO  TO  1001  002270 

81  IF(NNBR( I )-l)  83.83.65  002280 

83  CALL  STATEB(SVIP.VV.N.YY.L.8RC.D1.D2.PP.PM,SS.Q1.Q2.Q3.FINK)  002290 

65  CALL  FINBRA<U»VV.V.MAT.A,Z< I ) )  002300 
DO  66  J  =  1.6  002310 
DO  66  K  =  1.3  002320 

66  VV< J.K)=V< J.K)  002330 
IF(NNBR( I ) )68. 67,67  002340 

67  IF(NNBR(I ) -3 ) 1001 .68 »68  002350 

68  CALL  BRANCH(R3,VV,PHI ( I  I ) »U)  002360 
N  =  N+1  002370 
11=11+1  002380 
MAT-0  002390 
BC  =  1.  002400 
GO  TO  63  002410 

1001  CONTINUE  002420 
GO  TO  2000  002430 

C  002440 

C  002450 

C      CONTROL   OUT  OF  PLANE  002460 

C  002470 

998  Pl=1.0  002480 

P2=1.0  002490 

D1  =  0.  002500 

D2  =  0.  002510 

01  =  1.0  002520 

02=1.0  002530 

03=1.0  002540 

B1=0.  002550 

B2=0.  002560 

B3=0.  002570 

FINK=0.  002580 

IF(SVOd)  )1013. 1023. 1013  002590 

1013  CALL  STAVEO(SVOd)  .l.SVI  )  002600 

GO  TO  1043  002610 

1023  DO  1033  J=1.6  002620 

1033  SVK  J.l)=SVOP(  J.l)  002630 

1043  IF(SVO(NMBR) 11053.1063.1053  002640 

1053  CALL  STAVEO(SVO(NMBR) .l.SVE)  002650 

GO  TO  1083  002660 

1063  DO  1073  J=1.6  002670 

1073  SVE< J.l)=SVOP( J.NMBR)  002680 

1083  JK=1  002690 

IF(NBR)1002. 1002.1093  002700 

1093  DO  1095  N=2.NNN  002710 
IF(SVO(N) )1095, 1095. 1094  002720 

1094  CALL  STAVEO(SVOIN) .N.SVOP)  002730 

1095  CONTINUE  002740 

1002  N=2  002750 
11=1  002760 
HH=0.  002770 
DO  1003  1=1. NS  002780 
BC=0.  002790 


86 


MAT=0  002800 

P(  I  )=CTX( I )+CLZ( I )+CTY(  I  )  002810 

IF<P( I ) ) 1 .2.1  002820 

1  CALL  HANGEO(CTXU  >.CLZ(  I  )  »CTY(  I  >.U1  002830 
GO  TO  14  002840 

2  IF(AMYK)3.11.3  002850 

3  IF(RHO( I) )4,8.4  .  002860 

4  IF<  Z ( I ) )6.5.6  002870 

5  CALL  STIFOO  < THETA ( I > »RHO< I ) »U)  002880 
GO  TO  14  002890 

60CALL  CFIELO(AJT( I ) »R ( I ) »Z < I ) .G(  I ).SD.RHO< I ).THETA( I ) ,D< I > . DI ( I )  .FF   002900 

1AC.A)  002910 

7  CALL  POINO  <AMU(I).AIX( I ) .A  I Y ( I ) .OMEGA ( L > . TL ( I ) .Z U > ,RI • B J  002920 
CALL  MATMUL  (A.B.U)  002930 
MAT  =  1  002940 
GO  TO  14  002950 

8  IF(Z(I  )  J10.9.10  002960 

9  CALL  RIGIO  < AMU < I > .TL < I ) »A IX« I ) »OMEGA< L ) .AI Y < I > . AJT < I ) »RI »U )  002970 
GO  TO  14  002980 

100CALL  DISTMO(AMU(  I).AIXm.AIY(I)  .OMEGA  (L)  »TL<  I  )»AJT(  I ) »R< I ) »D< I > »D   002990 

II ( I ).G< I ) .SD.RI .FFAC.U)  003000 

GO  TO  14  003010 

11  IF(RHO( I) J4.12.4  003020 

12  IF(ZU)  U3.9.13  003030 

13  CALL  SFIELO  ( DI ( I ) .TL ( I ) . AJT ( I ) »R ( I ) .Z ( I ) »G<  I ) »SD.0( I  1 .FFAC  .  A)  003040 
GO  TO  7  003050 

14  IF(NNBRU))  201.200.201  003060 

200  IF(I-l)  15,202.15  003070 
202  CALL  STATEM(SVI.UU.D1.D2.P1.P2.PP.PM,SS)  003080 

GO  TO  15  003090 

15  CALL  FINMAT(U.UU.V.MAT.A»Z< I ) )  003100 
DO  16  J  =  1.6  003110 
DO  16  K=l,3  003120 

16  UU( J.K)»V< J.K )  003130 
IF(PM)  207.1114.85  003140 

1114  IF(FINK)  207.85.207  003150 

85  IF(BC)   207,207.1112  003160 

1112  IF(BRC)  207,207.208  003170 

208  CALL  BRCOR(R3.UU.XX,JK,VV,PM)  003180 
MMM=N-1  003190 
DO  209  J=1.3  003200 
DO  209  <=1,3  003210 

209  YY( J,K,MMM)=XX{ J.K)  003220 
207  GO  TO  1003  003230 

201  IF(NNBR( I )-l)   204,204,17  003240 
204  CALL  STATEB(SV0P,VV,N,YY,L,BRC,D1,D2,PP,PM,SS.Q1.Q2»Q3.FIN<)  003250 

17  CALL  FINBRA(U.VV»V.MAT.A.Z( I ) )  003260 
DO  18  J=1.6  003270 
DO  18  K=1.3  003280 

18  VV( J.K)=V( J.K)  003290 
IF(NNBR(I) 182.19.19  003300 

19  IF(NNBRU  )-3>1003,82.82  003310 
82  CALL  BRANCO(R3.VV.PHI( II ) »U)  003320 

N=N+1  003330 

11=11+1  003340 

MAT=0  003350 


37 


BC=  1. 

003360 

GO  TO  15 

003370 

1003 

CONTINUF 

003380 

c 

003390 

C 

ITERATION 

003400 

c 

003410 

2000 

CALL  DELMA  <  UU. SVE  »  V.  X2  t  Y2  t  KICK,  X  1.  FF,  PP.  HH.RRtSStYl.  REM,  F  I  NKtPM.Dl, 

003420 

LD2) 

003430 

IF(FINK)  87,87,86 

003440 

87 

IF(PM)  99,99,98 

003450 

98 

iF(KKK)  7030,7030,301 

003460 

301 

GO  TO  6001 

003470 

7030 

IF(FF-l.O)  7051,7060,7051 

003480 

7060 

KK=1 

003490 

M=M-1 

003500 

GO  TO  3000 

003510 

7051 

IF(HH-l.O)  300,6001,300 

003520 

300 

IF(SS)  303,302,302 

003530 

303 

DEL(L)=X2-Y2, 

003540 

D1=X1 

003550 

D2=(X2+Y2)/2.0 

003560 

GO  TO  61 

003570 

302 

DEL(L)=X1-Y1 

003580 

Dl=(Xl+Yl)/2.0 

003590 

D2  =  X2 

003600 

GO  TO  61 

003610 

990DEL(L)=V( 1 ,  1  ) »V < 2 . 2 ) *V < 3 ,3  -V< 1 , 3 > *V < 2 ,2 ) *V <  3, 1  +V( 1 .2 ) *V ( 2 t 3 >*V< 

003620 

,  1 )-V( 1,2)*V(2,1)*V(3,3)+V(1,3)*V( 3,2)*V(2,1  )-V( 1,1 )*V( 3,2)*V(2»3) 

003630 

GO  TO  61 

003640 

21 

PDEL=DEL(L) 

003650 

GO  TO  25 

003660 

22 

IF(AA)23,23,24 

003670 

23 

PDEL»0EL(L-1) 

003680 

GO  TO  2  5 

003690 

24 

PDEL=DELC 

003700 

25 

IF(PDEL*DEL(L) ) 26, 26 ,27 

003710 

26 

IF (AA) 28,28, 57 

003720 

27 

BB  =  0 

003730 

IF(AA)35.35,52 

003740 

28 

AA=1.0 

003750 

BB=1.0 

003760 

DELC*DEL(L-1) 

003770 

OMEGAC=OMEGA(L-l) 

003780 

GO  TO  59 

003790 

290OMEGAM ( M ) -0ME6AC+ ( OMEGA ( L ) -OMEGAC ) *ABSF ( DELC ) / ( ABSF ( DELC ) ♦ 

003800 

lABSFtOEL(L)  )  ) 

003  810 

GO  TO  32 

003820 

30 

OMEGA<L+1)=OMEGA(L)-0,9*DOM 

003830 

DOM=0.5*DOM 

003840 

31 

L  =  L+1 

003850 

BL  =  L 

003860 

IF(HL-BL)  7053,7053,86 

003870 

7053 

KA=1 

003880 

GO  TO  3000 

003890 

86 

IF(JK)1000, 1000, 1002 

003900 

32 

BM=M 

003910 

83 


33 

34 

340 


IF (HM-BM )3000.3000.3  3 
M=M+1 

IFIM-3)  340.34.34 

A00M= ( OMEGAM ( M-l ) -OMEGAM! M-2 ) ) «0. 1 
AA=0.0 
BB  =  0 
CC  =  0 
GG=0.0 
DOM=ADOM 
DELOO 

OMEGA (L+l )=OMEGAM(M-l)+2.*CR 
OMEGAO=OMEGA(L+l) 
GO  TO  31 

I F ( OMEGA ( L) -OMEGAG ) 5 3  .  60  .  60 
OMEGA ( L+l ) =OMEGA ( L) -( OMEGA ( L) -OMEGAC )/ 2 .0 
IF  (  OMEGA  (  L)  -OMEGA  <L+1)-CR>  51.51.  31 
CC=1.0 
GO  TO  31 
DELC=0EL(L) 
OMEGAC=OMEGA(L) 
OMEGA(L+l )*OMEGA(L)+DOM 
GO  TO  31 

IF(DELC*DEL(L) 129.29,55 
5  50OMEGAM( M ) =OMEGA ( L)  +  ( OMEGA ( L-l ) -OMEGA ( L)  ) *ABSF ( DEL  <  L )  ) / ( ABSF ( DEL ( L ) 
1>+ABSF(DEL(L-1) )) 
GO  TO  32 

I F ( OMEGA ( L ) -OME GAO ) 2 1 . 2 1 . 2 2 
IFJBBI58.58.50 
BB=1.0 

IF (OMEGA(L) -OMEGA ( L-l )-CR)29.?9.30 
M=M-1 

GO  TO  3000 
IF(GG)70.70.72 
IF(DEL(L))77.71.77 
GG=1.0 

0MEGA(L+1 )=OMEGA(L)-CR 
GO  TO  31 

IF(GG-1.0)73.73.74 
GG=2.0 
HDEL=DEL(L) 

OMEGA ( L+l )=OMEGA(L)+2.*CR 
GO  TO  31 

TDEL*DEL(L)*HDEL 
IF(TDEL)75.76.76 
OMEGAM ( M ) »OMEGA ( L-2 ) 
GO  TO  33 
KtC=l 
M-M-l 

GO  TO  3000 
IF(CC)56.56.54 


35 
50 

51 

52 

53 

54 


56 
57 
58 
59 
60 

61 

70 
71 


72 
73 


74 


75 


76 


77 


:  OUTPUT 

3000  IF(NM)400. 400.403 

400  PRINT  401 

401  FORMAT  UHl//> 


003920 
003930 
003940 
003950 
003960 
003970 
003980 
003990 
004000 
004010 
004020 
004030 
004040 
004050 
004060 
004070 
004080 
004090 
004100 
004110 
004120 
004130 
004140 
004150 
004160 
004170 
004180 
004190 
004200 
004210 
004220 
004230 
004240 
004250 
004260 
004270 
004280 
004290 
004300 
004310 
004320 
004330 
004340 
004350 
004360 
004370 
004380 
004390 
004400 
004410 
004420 
004430 
004440 
004450 
004460 
004470 


*? 


PRINT  402  004480 

402  FORMAT  <15H  PROGRAM  V I  PIPE . 35X ,  15HY.  S.  KIM  NHA-3.38X,  13H   MARCH  004490 

1  1966//.10X.108H               MODAL  FREQUENCIES  OF  A  PLANAR  PIPING  004500 

2  SYSTEM  ARE  DETERMINED  BY  AN  ITERATIVE  PROCEDURE  USING  THE  /29H  ME  004510 

3THOD  OF  TRANSFER  MATR I CES. /// 1 19H  * #♦»»»*#*»♦###•#**##**#•*##  004520 

4»#***##«**«*#*##*# *####*#*««*##**«»*»»*«*«*»**************** ******  004  5  30 

5##*#*»#»#»###**##«*     //)  00*540 

403  IF(NM)4058. 4058. 404  004550 

404  PRINT  405  004560 

405  FORMAT  <lHl//>  004570 

4058  NM=1  004580 

4059  IF( JO4060.4060.4065  004590 

4060  PRINT  4070. JJ  004600 
40700FORMAT  ( 50X .8HPROBLEM  I  2 » //35X  .47H IN  PLANE  MODE  FREQUENCIES  (RADIA  004610 

INS  PER  SECOND)    //45X .4HM0DE , 15X . 10H  FREQUENCY  /)  004620 

GO  TO  4079  004630 

4065  PRINT  4075. JJ  004640 

40750FORMATI 50X.8HPROBLEM  1 2 »/ /33X  ,5 1HOUT  OF  PLANE  MODE  FREQUENCIES  (RA  004650 

1DIANS  PER  SECOND)       //45X .4HM0DE . 15X . 10H  FREQUENCY       /)  00*660 

4079  MH=M  00*670 
PRINT  408.(M»OMEGAM(M) .M-l .MH)  004680 

408  FORMAT  (  46X  ,  I  2. 17X  .F14.8  /)  {J04690 
IF(KA)  5000,5000.5010  004700 

5010  PRINT  5020  004710 

5020  FORMAT(10X.53HPROBLEM  TERMINATED  DUE  TO  OMEGA  STORAGE  LIMITATION.  004720 

1   /)  004730 

5000  IF(KK)4080. 4090. 4080  004740 

4080  PRINT  4085  004750 
40850FORMAT(10X.69HPROBLEM  TERMINATED  DUE  TO  SIGNIFICANT  FIGURE  LIMITAT  004760 

HON  OF  COMPUTER.      /)  00*770 

4090  IF(AMYK)412,409,412  004780 

409  PRINT  411  004790 
4110FORMAT  (16X.53H               A  LUMPED  MASS  APPROACH  WAS  EMPLOYED.  A  00*800 

1ND  )  004810 

GO  TO  414  004820 

412  PRINT  413  004830 

413  FORMAT  (24X.45HA  DISTRIBUTED  MASS  APPROACH  WAS  EMPLOYED  AND       )  004840 

414  IF(SD)415.415.420  004850 

415  IF(RI )416.416.418  004860 

416  PRINT  417  004870 
4170FORMAT  (2X.67HTHE  EFFECTS  OF  SHEAR  DEFLECTION  AND  ROTARY  INERTIA  W  004880 

1ERE  CONSIDERED.  /)  004890 

GO  TO  4610  004900 

418  PRINT  419  004910 

4190FORMAT  (3X.63HTHE  EFFECTS  OF  SHEAR  DEFLECTION  WERE  CONSIDERED  WHIL  004920 

IE  THOSE  OF      /3X.35HROTATIONAL  INERTIA  WERE  NEGLECTED.        /)  004930 

GO  TO  4610  004940 

420  IF ( RI )421. 421.423  004950 

421  PRINT  422  004960 
4220FORMAT  (3X.72HTHE  EFFECTS  OF  ROTATIONAL  INERTIA  WERE  CONSIDERED  WH  004970 

1ILE  THOSE  OF  SHEAR           /3X .27HDEFLECT ION  WERE  NEGLECTED.    /)  004980 

GO  TO  4610  004990 

423  PRINT  424  005000 

4240FORMATI3X.66HTHE  EFFECTS  OF  SHEAR  DEFLECTION  AND  ROTARY  INERTIA  WE  005010 

IRE  NEGLECTED.  /)  005020 

4610  IF(PM)  425.425.4713  005030 


90 


4713  IF(SS)  4611,4612.4612  005040 

4611  IF(PP-1.)  4613.4614,4615  005050 

4613  PRINT  4623  005060 

4623  FORMAT!24X,35HPESTEL  MAHRENHOLTZ  METHOD(A)  USED.     /)  005070 
GO  TO  425  005080 

4614  PRINT  4624  005090 

4624  FORMAT(24X.35HPESTEL  MAHRENHOLTZ  METHOD(B)  USED.     /)  005100 
GO  TO  425  005110 

4615  PRINT  4625  005120 

4625  FORMAT(24X.35HPESTEL  MAHRENHOLTZ  METHOD(C)  USED.     /)  005130 
GO  TO  425  005140 

4612  IF(SS-1.)  4616,4617,4618  005150 

4616  PRINT  4626  005160 

4626  FORMAT! 24X,35HPESTEL  MAHRENHOLTZ  METHOD(D)  USED.     /)  005170 
GO  TO  425  005180 

4617  PRINT  4627  005190 

4627  FORMAT! 24X.35HPESTEL  MAHRENHOLTZ  METHOO(E)  USED.     /)  005200 
GO  TO  425  005210 

4618  PRINT  4628  005220 

4628  FORMAT(24X,35HPESTEL  MAHRENHOLTZ  METHOD(F)  USED.     /)  005230 
425  IF(PM)  9010,9020,9000  005240 

9010  PRINT  9012  005250 

9012  FORMAT! 24X.20HDELTA  METHOD  USED.    /)  005260 

GO  TO  9000  005270 

9020  PRINT  9022  005280 

9022  FORMAT! 24X , 29HMODI FI ED  DELTA  MFTHOD  USED.     /)  005290 

9000  PRINT  440  005300 

4400FORMAT! H8H  ******************************************************   005310 

1 *****»******»**#**»*#*****»*#♦#**»******#*♦»**#***#**»«#«  **»##  )   005320 

IF(GDABC-1. 14511,4520,4401  005330 

4401  PRINT  4402  005340 

44020FORMATI35H  SECTION  LENGTH  AND  SUBSECTION  DATA/23X, 14HSECT ION  NUMBE   005350 

1R,8X,25HLENGTH  OF  SECT  ION ( INCHES ), 5X, 2 1HNUMBER  OF  SECTIONS  )    005360 

PRINT  443, ( I.TL! I ) »Z < I ).I»1,NS  )  005370 

443  FORMAT  ! 29X » 13, 23X ,F7.2 ,21X ,F4.0 )  005380 
LM=L  005390 
IF(PM)  7054,7054,7055  005400 

7055  PRINT  441  005410 

4410FORMAT  (11H  GRAPH  DATA  /22X , 16HI TERATION  NUMBER, 4X .30HFREQUENCY  (R   005420 

1ADIANS  PER  SECOND) »3X,30H    VALUE  OF  REMAINDER  )    005430 

GO  TO  705  7  005440 

7054  PRINT  7056  005450 

70560FORMAT  ( 1 1H  GRAPH  DATA  /22X, 16HI TERATION  NUMBER ,4X,30HFREOUENCY  (R   005460 

1ADIANS  PER  SECOND) ,3X,30HVALUE  OF  FREQUENCY  DETERMINANT  )    005470 

7057  PRINT  444, (L, OMEGA (L) .DEL < L ) ,L«1 »LM)  005480 

444  FORMAT  ( 29X , I 3,23X,F1 1  .6, 18X ,E15.8  )  005490 
4520  PRINT  442  005500 

442  FORMAT! IX, 10HINPUT  DATA  /)  005510 

PRINT  4521  005520 

45210FORMAT! IX ,9HCOMPONENT ,4X,8HDI AMETER, 5X , 14HWALL  THICKNESS  »5X,6HLEN   005530 

1GTH,5X,9HRADIUS  OF ,9X ,8H I NCLUDED,5X , 7HDENS ITY , 5X , 7HELAST IC ,5X ,5HSH   005  540 

2EAR  /57X,9HCURVATURE»5X.12HANGLE  OF  ARC. 16X, 7HMODULUS,5X .7HMODULUS   005550 

3  /)  005560 

DO  4525  1  =  1, NS  005570 

IF(P(I ) 14524,4524,4523  005580 

4524  PRINT  4522 , I ,D(  I )  ,DT ( I  )  ,TL ( I ) ,RHO( I ) ,THETR ! I ) ,BAMU( I ) ,E ( I ) ,G< I )  005590 
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4  5220FORMAT  (3X,I3»8X,F7.3»10X,F6.3»7X,F8.3»4X,F8.3.9X.F6.2.8X.F7.2.4X.   005600 

1E8.2.4X.E8.2)  005610 

GO  TO  4525  005620 

4523  SH=1.  005630 

4525  CONTINUE  005640 
IF(SH)4530. 4530. 4526  005650 

4526  PRINT  4532  005660 
4532  FORMAT  ( / /1X .7HHANGERS  /)  005670 

PRINT  4527  005680 

452  70FORMAT  (2X.9HCOMPONENT.14X.3HCLX.14X.3HCLY.14X.3HCLZ.14X.3HCTX.14X   005690 

1.3HCTY.14X.3HCT2   /)  005700 

4528  DO  4530  1*1. NS  005710 
IF(P(I ) )4530. 4530. 4529  005720 

4529  PRINT  4531 . I .CLX ( I ) . CLY ( I ) ,CLZ ( I ) »CTX ( I ) .CTY < I ) .CTZ ( I )  005730 
4  531  FORMAT(3X.I3H5X.F10.2»7X»F10.2.7X.F10.2.7X,F10.2.7X,F10.2.7X.F10.   005740 

12  )  005750 

4530  CONTINUE  005760 
PRINT  448  005770 

448  FORMAT  J//20H  BOUNDARY  CONDITIONS  )  005780 
PRINT  449,  (SVK  1,1  ).  1*1. 6)  005790 

449  FORMAT  (5H  SVI=6F3.0>  005800 
PRINT  451,(SVE( 1.1  ). 1  =  1.6)  005810 

451  FORMAT  (5H  SVE=6F3.0)  005820 

IF(NBR)4511. 4511, 4502  005830 

4502  IF(JK)4505. 4505. 4515  005840 

4505  PRINT  4510, ( (SVIP( I, J) ,1=1,6) ,J=2,NNN)  005850 

4510  FORMAT  (6F3.0)  005860 
GO  TO  4511  005870 

4515  PRINT  4517,( (SVOPI  I.J)  .1  =  1.6) .J*2.NNN)  005880 

4517  FORMAT  (6F3.0)  005890 

4511  CONTINUE  005900 
C      SORT  005910 

IFJGRAPH-1.)  6001.6000.6001  005920 

6000  LM1*LM-1  005930 

DO  6100  L*1.LM  005940 

EDL(L)-DEL(L)  005950 

6100  CONTINUE  005960 

DO  6600  L=1.LM1  005970 

LP1*L+1  005980 

DO  6600  M*LP1,LM  005990 

IF (OMEGA(L) -OMEGA ( M) )  6600,6600,6003  006000 

6003  TEMP*OMEGA(L)  006010 

OMEGAtL)"OMEGA(M>  006020 

OMEGA(M)=TEMP  006030 

TEMP=EDL(L)  006040 

EDL(L)*EDL(M)  006050 

EDL(M)=TEMP  006060 

6600  CONTINUE  006070 

C      GRAPH  006080 

NUMPTS=60  006090 

<=0  006100 

DO  6006  J=l,12  006110 

6006  ITITLE(J)=8H  006120 

ITITLE(3)=8H  KIM  Y  S  006130 

ITITLE(8)=8HFREQUENC  006140 

ITITLE(9)=8HY   VS   R  006150 


n 


6020 
6021 


ITITLEl 10)=8HEMAINOER 
LT1=LM/NUMPT5 
IF(LTl)  6024,6025.6024 
6025  NUMPTS=LM 
LABEL=4H 

cZlL  DRAWINUMPTS, OMEGA, EDL, 0,0, LABEL, ITITLE, 0,0, 0.0. 0.0. 9. 

110. K. LAST) 
GO  TO  6001 
6024  DO  6020  J=l»17  ■ 
LT2  =  0 

LT2=LT2+J*NUMPTS 
IF(LT2-LM)  6020.6021.6022 
CONTINUE 
LT=LT1 
GO  TO  6023 

6022  LT=LT1+1 
NUMEND=LM-LT1*NUMPTS 

6023  L=l 

DO  7021  1  =  1. LT 
IF(I-U  6008.6041.6008 
IF(I-LT)  8011.6011.6011 
IF(I-3)  6042.6043,8012 
IP ( i-5 )  6044,6045,8013 
IF(I-7)  6046,6047,8014 
I F ( i-9 )  6048,6049,8015 
IF(I-LT)  6050,6011,6001 
LABEL=4H  ONE 
GO  TO  6009 
LABEL=4H  TWO 
GO  TO  6010 
LABEL=4HTHRE 
GO  TO  6010 
LABEL=4HF0UR  . 
GO  TO  6010 

6045  LABEL*4HFIVE 
GO  TO  6010 

6046  LABEL=4H  SIX 
GO  TO  6010 
LABEL=4HSEVN 
GO  TO  6010 
LABEL=4H  EIT 
GO  TO  6010 
LABEL=4HNINE 
GO  TO  6010 
LABEL=4H  TEN 
GO  TO  6010 

6009  DO  6031  N*1.NUMPTS 

EDL«N)=EDL(L+N-1) 
6031  OMEGA(N)*OMEGA<L+N-l) 

FO^ATUHW^^^  .  o 

CALL  DRaSInUMPTS^OMEGA.EDL.O.O.LABEL. IT. TLE.0,0.0.0. 0.0.9. 

110.K.LAST) 

GO  TO  7021 

6010  L=L+NUMPTS 


6008 
8011 
8012 
8013 
8014 
8015 
6041 

6042 

6043 

6044 


6047 
6048 
6049 
6050 


6027 


006160 

006170 

006180 

006190 

006200 

006.210 

006220 

006230 

006240 

006250 

006260 

006270 

006280 

006290 

006300 

006310 

006320 

006330 

006340 

006350 

006360 

006370 

006380 

006390 

006400 

006410 

006420 

006430 

006440 

006450 

006460 

006470 

006480 

006490 

006500 

006510 

006520 

006530 

006540 

006550 

006560 

006570 

006580 

006590 

006600 

006610 

006620 

006630 

006640 

006650 

006660 

006670 

006680 

006690 

006700 

006710 
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OMEGAX=OMEGA(L)  006720 

DO  6030  N=1.NUMPTS  006730 

EDL(N)=EDL(N+L-1)  006740 

OMEGA(N)=OMEGA(L+N-l )  006750 

6030  OMEGA(N)=OMEGA(N)-OMEGAX  006760 

PRINT  6032. NUMPTS.EDU  1  1  006770 

6032  FORMAT(///////5X,7HNUMPTS=.  13   .  5X . 7HEDL ( 1 ) = . E 12  .9  )  006780 
CALL  DRAW(NUMPTS.OMEGA.EDL»0.0.LABEL. ITITLE. 0.0.0. 0.0. 0.9.  006790 

110.K.LAST)  006800 

GO  TO  7021  006810 

6011  L=L+NUMPTS  006820 

OMEGAX*OMEGA(L>  006830 

DO  6033  N*1,NUMEND  006840 

EDL(N)=EDL(N+L-1)  006850 

0MEGA«N)=0MEGA(N+L-1>  006860 

6033  0MEGA(N)=0ME5A(N)-0MEGAX  006870 
PRINT  6051.NUMEND,EDL(1)  006880 

6051  FORMAT(///////5X,7HNUMEND=.I3.    5X  .7HEDL( 1 ) * . E12 . 9  )  006890 

LABEL»4HLAST  006900 

CALL  DR AW (NUMEND. OMEGA. EDL. 0.0. LABEL, ITITLE. 0.0. 0.0. 0.0. 9.  006910 

110. K, LAST)  006920 

GO  TO  6001  006930 

7021  CONTINUE  006940 

6001  IF(HH-l.O)  8024,8023.8023  006950 

8023  MH=M  006960 
LM«L  006970 

8024  DO  446  M=1.MH  006980 
446  OMEGAM(M)*0.0                                                        006990 

DO  447  L*1»LM  007000 

OMEGA(L)=0«0  007010 

EDL(L)«0.  007020 

007030 
007040 
007050 
007060 
007070 
007080 
007090 
007100 
007110 
007120 
007130 
007140 
007150 
007160 
007170 
007180 
007190 
007200 
007210 
007220 
007230 
007240 
'  007250 

007260 
007270 


447 

DEL(L)»0.0 

DELC=0.0 

OMEGAC-0.0 

IF(HH-l.O)  8001,8022,8001 

8001 

IF(PM)  8031,8031,8033 

8033 

IF(REM)  8031,8032,8032 

8031 

IF(MN)8022, 8022, 452 

8032 

IF(REM-1.)  8041,8042,8043 

8041 

IF(PP-l)  8043,8043,8072 

8072 

IF(SS)  8073,8074,8074 

8073 

ss«o. 

GO  TO  8022 

8074 

IF (SS-1. 18042, 8042, 8075 

8075 

SS=-1. 

PP-O. 

REM=-1, 

GO  TO  8033 

8042 

SS=SS-1. 

IF(SS)  8051,8052.8053 

8051 

SS=1. 

GO  TO  8022 

8052 

SS  =  2. 

GO  TO  8022 

8053 

SS=-1. 

REM=-1. 

% 


GO  TO  8033 

8043 

pp=pp-l. 

IF(PP)  8061»806?.8063 

8061 

PP  =  1. 

GO  TO  8022 

8062 

PP  =  2. 

GO  TO  8022 

8063 

PP  =  0. 
REM=-1. 
GO  TO  8033 

452 

PP  =  0. 
RR*0. 
SS=-1. 

8022 

L»l 

M=l 

AA=0.0 

BB=0.0 

CC=0.0 

FF=0.0 

GG=0.0 

KK  =  0 

KA  =  0 

KKK  =  0 

ADOM=BDOM 

CR  =  BR 

OMEGA! 1 )=OMEGAI 

OMEGAO»OMEGAI 

DOM=ADOM 

IFIHH-1.0)  8021.7040.8021 

7040 

IF(SS)  8191.8192.8192 

8191 

IF(PP-1.)  8091.8092*8093 

8192 

IF(SS-1.)  8291.8292.8293 

8091 

PRINT  8081 

8081 

F0RMAT(1H1.////18HMETH0D   A 
GO  TO  8021 

FAILED. 

8092 

PRINT  8082 

8082 

FORMAT! 1H1.////18HMETHOD   8 
GO  TO  8021 

FAILED. 

8093 

PRINT  8084 

8084 

FORMATt 1H1.////18HMETHOD   C 
GO  TO  8021 

FAILED. 

8291 

PRINT  8281 

8281 

FORMAT! 1H1.////18HMETHOD   D 
GO  TO  8021 

FAILED. 

8292 

PRINT  8282 

8282 

FORMAT! 1H1.////18HMETHOD   E 
GO  TO  8021 

FAILED. 

8293 

PRINT  8283 

8283 

FORMAT! 1H1.////18HMETHOD   F 

FAILED. 

8021 

IF(PM)  453.453.8035 

8035 

IF(REM)  8037.8036.8036 

8036 

IF(JK)  998.999.998 

8037 

REM=ERM 

453 

IF!PM)  6081.6082.6080 

6081 

IFIPMM-1.)  6084.6080.6084 

6084 

IF!PMM*0. 1-10.1)  6086.6085.6086 

///) 


///) 


///) 


///) 


///) 


///) 


007280 
007290 
007300 
007310 
007320 
007330 
007340 
007350 
007360 
007370 
007380 
007390 
007400 
007410 
007420 
007430 
007440 
007450 
007460 
007470 
007480 
007490 
007500 
007510 
007520 
007530 
007540 
007550 
007560 
007570 
007580 
007590 
007600 
007610 
007620 
007630 
007640 
007650 
007660 
007670 
007680 
007690 
007700 
007710 
007720 
007730 
007740 
007750 
007760 
007770 
007780 
007790 
007800 
007810 
007820 
007830 
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007840 

6085  PM-.1.  00lll°n 
GO  TO  8036  007860 

6086  Pm=0.  007870 
GO  TO  8036  •  ;                                        007880 

6082  IF(PMM*0. 1-1.1)  6080,6080.6085  007890 

6080  !F(MN)  9001 .9001 .8034  007900 

8034  MN  =  0  007910 

PM=PPM  007920 

GO  TO  998  007930 

9001  JJ=JJ+1  007940 

DO  454  I=1.NS  007950 

P(  I  )"0.  007960 

454  Z( I  )=0.0  007970 

GO  TO  450  007980 

600  END  007990 


C 


008000 


SUBROUTINE  SUBSEC  ( RHO .THETA ,D. TL .HM.Z .OMEGA  I ,HMH)  008010 

RHOD=ABSF(RHO)  008030 

IF(HMH)  41.40.41  008040 

40  IFIHM-5. 0)8.8*7  008050 
8  IF(OMEGAI-800.)  10.7.7  008060 
7  HM=6.  008070 

GO  TO  10  008080 

41  HM=HMH  008090 
10  RA=TL/0  008100 

IF(RHOD-0. )3.3.4  008110 

3  IF(RA-1.)5.5.6  008120 

5  Z=0.  008130 
RFTURN  008140 

6  IF(RA-3.)  15.15.16  008150 

15  Z=2.0  008160 
RETURN  008170 

16  IF(RA-6.) 17.17.18  008180 

17  Z«3.0  008190 
RETURN  008200 

18  IFCHM-3.) 1.1.2  008210 

1  Z«6.0  008220 
RETURN  008230 

2  Z=2.0*HM  008240 
RETURN  008250 

4  RD=RHOD/D  008260 

IF(RD-1.)19.19,23  008270 

19  Z=0  008280 

RETURN  008290 

23  IF(RD-3. 124.24.25  008300 

24  IF(THETA-.44)26.26.27  008310 

26  Z=0.  008320 
RETURN  008330 

27  IF(THETA-1.57)28.28.29  008340 

28  Z=4.0  008350 
RETURN  008360 

29  Z=6.0  008370 
RETURN  008380 

25  IF(RD-6. 130.30. 31  008390 

30  IF(THETA-.44)32,32.33 
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32  Z  =  4.0 

RETURN 

33  IF(THETA-1.57)34,34,35 

34  Z  =  6.0 
RFTURN 

35  Z=8.0 
RFTURN 

31  IF(THETA-.44)36.36.3 

36  Z  =  4.0 
RETURN 
END 


SUBROUTINE  DISTM  ( AMU. A  I Y.OMEGA.TL .R.D.DI »G.SD»RI .E t FFAC.U ) 

DIMENSION  U(6»6) 

ARA=( 3. 141592654/4.0 ) * ( D**2-DI »*2 > 

AREA=ARA«FFAC 

BE=TL*0MEGA*SQRTF(AMU/(E*ARA1 ) 

PA=R*AMU»TL**4*OMEGA**2 

IFIRI )1.1  .2 

1  T=R*AMU*(AIY*0MEGA*TL)**2 
GO  TO  3 

2  T=0.0 

3  IF<SD)4.4,5 

4  S=(AMU*(OMEGA*TL)**2)/(G*AREA) 
GO  TO  6 

5  5=0.0 

6  SPT=S+T 
VV=SQRTF(P4+(S-T)**2/4.0) 
AL1=SQRTF(VV-.5*SPT) 
AL2=SQRTF(VV+.5»SPT) 
CL=1.0/(AL1**2+AL2*»2) 
CHL1=(EXPF( AL1)+1./EXPF(AL1) )/2. 
SHL1  =  (EXPF{AL1)-1./EXPF(ALU  )/?. 
CL2=COSF(AL2> 

SL2=SINF(AL2) 

SBE=SINF(BE)  " 

CBE=COSF(BE) 

COO=CL*(AL2**2*CHL1+AL1**2*CL2) 

C01=CL*(  ( AL2**2*SHL1  )  /  AL  1  +  ( AI_1#*2*5L2  )  /  AL2  ) 

C02=CL*(CHL1-CL2) 

C03=CL*«SHL1/AL1-SL2/AL2) 

DO  7  J=1.6 

DO  7  K=1.6 

7  U(J 
U(l 
U(  1 
U(2 
U<2 
U(2 
U(2 
U(3 
U(3 
U(3 
U(3 
RR  = 


K)=o.n 

1)=»CBE 

6)=TL*SBE/(BE*ARA»E) 
2)=CO0-S*CO2 
3)=TL*(C01-SPT*C03) 
4)=R»C02*TL«*2 

5)=(R*TL»«3/P4)»( (P4+S**2l*C03-S*C01) 
2)=P4*C03/TL 
3)=COO-T*C02 
4)=TL*R«(C01-T*C03) 
5)=U(2.4) 
l.O/R 


008400 
008410 
008420 
008430 
008440 
008450 
008460 
008470 
008480 
008490 
008500 
008510 
008520 
008530 
008540 
008550 
008560 
008570 
008580 
008590 
008600 
008610 
008620 
008630 
008640 
008650 
008660 
008670 
008680 
008690 
008700 
008710 
008720 
008730 
008740 
008750 
008760 
008770 
008780 
008790 
008800 
008810 
008820 
008830 
008840 
008850 
008860 
008870 
008880 
008890 
008900 
008910 
008920 
008930 
008940 
008950 
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U(4,2)=P4*RR*C02/TL**2 

U(4.3)=RR*< (P4+T*»2)*C03-T#C01 ) /TL 

U(4.4)=U( 3»3) 

U(4»5)=U(2»3) 

U(5.2)=P4*RR*(C01-S*C03)/TL«*3 

U(5»3)=U<4,2) 

U(5»4)=U<3.2) 

U(5,5)=U<2.2) 

U(6»t)=-AMU*TL*OMEGA**2*SBE/BE 

U(6»6)=CBE 

RETURN 

END 

c 

c 

SUBROUTINE  DISTMO  < AMU  .A  I  X tA I Y .OMEGA . TL . A JT .R .D.D I .G.SD.RI .FFAC.U) 

DIMENSION  U16.6) 

W=1./(G*AJT) 

AREA=( 3.141592654/4.0 )*<D*»2-DI»*2)*FFAC 

BE=+SORTF(AMU  *  (  T  L   *AI X*OMEGA ) **2*W ) 

P4=R*AMU    »  TL  »«4*OMEGA**2 

IF(RI)2.2.3 

2  T=R«AMU»(AIY»OMEGA*TL)»«2 
GO  TO  4 

3  T=0.0 

4  IF(SD)5.5.6 

5  S=(AMU*(OMEGA*TL)**2 )/ (G»AREA) 

GO  TO  7 

6  S  =  0.0 

7  VV=+SQRTF(P4+(S-T)*»2/4.0) 
SPT=S+T 

AL1=+SQRTF(VV-.5»SPT) 
AL2=+SQRTF(VV+.5»SPT> 
CL*1.0/(AL1**2+AL2**2) 
CHL1*<EXPF(AL1)+1./EXPF(AL1) )/2. 
SHL1*(EXPF(AL1)-1./EXPF(AL1> )/2. 
CL2»COSF(AL2) 

SL2«SINF(AL2)  - 
SBE*SINF(BE)  ■ 
CBE=COSF(BE> 

CO0*CL*(AL2**2*CHLl  +  AL1**2*CL2) 
C01«CL*( <AL2»*2*SHL1 1/AL1  +  ( ALl**2*SL2 ) /AL2 ) 
C02=CL*(CHL1-CL2> 
C03=CL*(SHL1/AL1-SL2/AL2) 
DO  1  J«1.6 
DO  1  K«1.6 
1  U(J.K>«0.0 
U(1.1)=CBE 
U(1.2)=TL*W*SBE/BE 

U(2.1)=-AMU   *  TL   *AIX**2«SBE/BE 
U(2,2)»CBE 
U(3.3)=CO0-S*CO2 
U<3.4)=  TL  *(C01-SPT*C03) 
U(3.5)=R»C02*TL*»2 

U(3.6)=(R»TL**3/P4)*(-S*C01+(P4+S»*2)*C03) 
U(4,3)=P4*C03/TL 
U(4,4)=CO0-T»CO2 


008960 
008970 
008980 
008990 
009000 
009010 
009020 
009030 
009040 
009050 
009  060 
009070 
009080 
009090 
009  100 
009110 
009120 
009130 
009  140 
009150 
009160 
009170 
009180 
009190 
009200 
009210 
009220 
009230 
009240 
009250 
009260 
009270 
009280 
009290 
009300 
009310 
009320 
009330 
009340 
009350 
009360 
009370 
009380 
009390 
009400 
009410 
009420 
009430 
009440 
009450 
009460 
009470 
009480 
009490 
009500 
009510 


n 


30 


U(4,5)=  TL*R*(C01-T«C03) 

U(4.61=U<3.5) 

RR=1.0/R 

UC5.3)=P4*RR*C02/TL*#2 

U(5.4)=RR*(-T*C01+(P4+T**2)*C03)/TL 

U(5.5)=U<4,4) 

U(5.6)=U(3.4) 

U(6.3)=P4*RR»(C01-S*C03)/TL**3 

U(6.4)=U( 5.3) 

U(6.5)=U<4,3) 

U(6.6)=U(3.3) 

RETURN 

END 


SUBROUTINE  RIGID  ( AMU »TL. OMEGA »A I Y ,R I »U ) 

DIMENSION  U<6.6) 

AM=AMU*TL 

DO  1  J=1.6 

DO  1  K  =  l»6 

U(J.K)=0.0 

U(l.l)=1.0 

U(6.1)=-AM*0MEGA**2 

U(6.6)=1.0 

U(2.2)=1.0 

U(2.3)=TL 

U(3»3)-1.0 

U(4.2)=AM»TL*OMEGA**2/2.0 

IF(RI)2.2»3 

U(4,3)=AM*0MEGA**2*< TL»*2/6.0-A I Y**2 ) 

GO  TO  4 

U(4»3)=AM*(OMEGA*TL)»*2/6.0 

U(4.4)*1.0 

U(4,5)=TL 

U(5,2)=-U(6.1) 

U(5.3)=U<4,2> 

U(5,5)=1.0 

RETURN 

END 

SUBROUTINE  RIGIO  ( AMU.TL. A  IX .OMEGA . AI Y .AJTtSDRI .U ) 

DIMENSION  UI6.6) 

DO  2  J=1.6 

DO  2  K=l»6 

U(J.K)=0.0 

AM=AMU*TL 

Ud.lW.O 

U(  2  » 1 ) *-AMU*TL* ( A  I X*OMEGA ) **2 

U(2.2)=1.0 

U(3.3>=1.0 

U(3.4)=TL 

U(4.4)*1.0 

U(5.3)=AM*TL*OMEGA«*2/2.0 

IF(SDRI-0.)  30.30.31 

U(5.4)=AM*0MEGA**2*<  TL**2/6. 0-A I Y**2 > 


009520 
009530 
009540 
0095  50 
009  560 
009570 
009580 
009590 
009600 
009610 
009620 
009630 
009640 
009650 
009660 
009670 
009680 
009690 
009700 
009710 
009720 
009730 
009740 
009750 
009760 
009770 
009780 
009790 
009800 
009810 
009820 
009830 
009840 
009850 
009860 
009870 
009880 
009890 
009900 
009910 
009920 
009930 
00994G 
009950 
009960 
009970 
009980 
009990 
010000 
010010 
010020 
010030 
010040 
010050 
010060 
010070 


99 


GO  TO  1 
31  U<5.4) =AM*OMEGA**2»( Tl*»2/6.0) 
1  U(5.5)  =  ]  .0 
U(5,6)=TL 

U(6.3)=AM*0MFGA»*2 
U<6.4)=U<5.3) 
U(6»6)=1.0 
RFTURN 
END 


SUBROUTINE  STIFCO  ( THETA tRHO tU ) 
DIMENSION  UI6.6) 
DO  1  J=l»6 
DO  1   K*l»6 

1  U(J.  10=0.0 
!F(RHO)2»2.3 

2  V=-1.0 
GO  TO  4 

3  V=1.0 

4  CT=COSF(THETA) 
ST=SINF(THETA) 
U(1#1)«CT 
U(1.2)»-V*ST 
U(2.1)»V*ST 
U(2.2)=CT 
U(3.3)»1.0 
U<4,4)=1.0 
U(5.5)«CT 
U(5.6)=-V*ST 
U(6.5)=V*ST 
U(6.6)*CT 
RETURN 

END 


SUBROUTINE  STIFOO  < THETA .RH0»U ) 

DIMENSION  U(6t6) 

IF(RHO) 1*1*2 

1 

V=-1.0 

GO  TO  3 

2 

V=1.0 

3 

CT=COSF(THETA) 

ST=SINF(THETA) 

DO  4  J=l»6 

DO  4  K  =  1.6 

4 

U( J.K)=0.0 

U(1.1)=CT 

U(  1.4)=ST*V 

U(2.2)=CT 

U(2.5)=ST*V 

U(3,3)=1.0 

U(4»1)»-ST*V 

U(4,4)=CT 

U(5»2)=-ST*V 

U(5,5)=CT 

010080 

010090 

010100 

010110 

010120 

010130 

010140 

010150 

010160 

010170 

010180 

010190 

010200 

010210 

010220 

010230 

010240 

010250 

010260 

010270 

010280 

010290 

010300 

010310 

010320 

010330 

010340 

010350 

010360 

010370 

010380 

010390 

010400 

010410 

010420 

010430 

010440 

010450 

010460 

010470 

010480 

010490 

010500 

010510 

010520 

010530 

010540 

010550 

010560 

010570 

010580 

010590 

010600 

010610 

010620 

010630 


I 


010640 
010650 
010660 
010670 
010680 
010690 
010  700 
010710 
010720 
010730 
010740 
010750 
.  010760 
010770 
010780 
010790 
010800 
010810 
010820 
010830 
010840 
010850 
010860 
010870 
010880 
010890 
010900 
010910 
010920 
010930 
010940 
010950 
010960 
010970 
010980 
010990 

SUBROUTINE  STAVEO  (SV.N.BC)  011000 

DIMENSION  BC(6.22)  011010 

DO  10  L=l»6  011020 

10  BC(LtN)=0.0  011030 

IF(SV-1.0)1,1.2  011040 

1  BC(2.N)=1.0  011050 
BC(5.N)=1.0  011060 
BC(6.N)=1.0  011070 
RETURN  011080 

2  IF(SV-2.0)3,3.4  011090 

3  BC(ltN)=1.0  011100 
BC(3,N)=1.0  011110 
BC(4tN)=1.0  011120 
RETURN  011130 

4  IF<SV-3.0)5,5.6  011140 

5  BC(2.N)=1.0  011150 
BC(4,N)=1.0  011160 
BC(6.N)=1.0  011170 
RETURN                                                               011180 

6  IF(SV-4.0)7f7.5  011190 


U(6t6)=1.0 

RETURN 

END 

c 

c 

SUBROUTINE  ST AVEC < SV .N .BC ) 

DIMFNSION  BC<6,22> 

DO  10  J=l,6 

10 

BC( J.N)=0. 
IF(SV-1.)1.1.2 

1 

BC(4.N)=1.0 
BC(5.N)=1.0 
BC(6tN)=1.0 
RETURN 

2 

IF(SV-2.)3»3.4 

3 

BC(1.N)=1.0 
BC(2»N)=1.0 
BC(3.N)=1.0 
RETURN 

4 

IF(SV-3.)5»5,6 

5 

BC(3»N)=1.0 
BC(5tN)=1.0 
BC(6»N)=1.0 
RETURN 

6 

IF(SV-4.)7»7.8 

7 

8C(liN)=1.0 
BC(3»N)=1.0 
BC(5.N)=1.0 
RETURN 

8 

BC(2.N)=1.0 

BC(3.N)=1.0 

BC(6.N)=1.0 

RETURN 

END 

c 

c 

7 

BC(1.N)=1.0 

BC(4.N)=1.0 

BC(6.N)=1.0 

RETURN 

END 

c 

c 

SUBROUTINE  INVERT  (A.N.D.L.M) 

c 

PROGRAM  FOR  FINDING  THE  INVERS 

c 

DIMENSION  AI25.25) »L( 25) .M<25) 

c 

SEARCH  FOR  LARGEST  ELEMENT 

0=1.0 

DO80  K=1.N 

L(K)=K 

M(K)=K 

BIGA  =  A(K..K) 

DO20  I=K.N 

DO20  J  =  K.N 

IF(ABSF(BIGA)-ABSF(A{ I ,J) ) )  10 

10 

BIGA=A( I.J) 

L(K)=I 

M(K)=J 

20 

CONTINUE 

c 

INTERCHANGE  ROWS 

J=L(K) 

IF(LIKJ-K)  35.35.25 

25 

DO30  1=1. N 
HOLD=-A(K.I ) 
A(K,I )»A( J.I ) . 

30 

A(  J.I )=HOLD 

c 

INTERCHANGE  COLUMNS 

35 

I=M(K) 

IF(M(<)-K)  45.45.37 

37 

D040  J=1.N 
HOLD=-A(J,K) 
A< J.K)*A( J.I  ) 

40 

A( J.I )=HOLD 

c 

DIVIDE  COLUMN  BY  MINUS  PIVOT 

45 

D055  I-1»N 

46 

IF( I-K)50.55.50 

50 

A(  I  »IC)=A{  I.K)/(-A(K.K)  ) 

55 

CONTINUE 

c 

REDUCE  MATRIX 
D065  1=1. N 
D065  J=1.N 

56 

IF(I-K)  57,65.57 

57 

IF(J-IC)  60.65.60 

60 

A(I.J)=A(I.K)*A(IC.J)+A{I.J) 

65 

CONTINUE 

c 

DIVIDE  ROW  BY  PIVOT 
D075  J=1.N 

68 

IF( J-K)70.75.70 

70 

A(K.J)=A(K»J)/A(K.K) 

75 

CONTINUE 

c 

CONTINUED  PRODUCT  OF  PIVOTS 

011200 
011210 
011220 
011230 
011240 
011250 
011260 
011270 
011280 
011290 
011300 
011310 
011320 
011330 
011340 
011350 
011360 
011370 
011380 
,20.20  011390 

011400 
011410 
011420 
011430 
011440 
011450 
011460 
011470 
011480 
011490 
011500 
011510 
011520 
011530 
011540 
011550 
011560 
011570 
011580 
011590 
011600 
011610 
011620 
011630 
011640 
011650 
011660 
011670 
011680 
011690 
011700 
011710 
011720 
011730 
011740 
011750 


02 


D=D*A(K,K)  011760 

C      REPLACE  PIVOT  BY  RECIPROCAL  011770 

A(K,K)=1.0/A(K,KJ  011780 

80  CONTINUE  011.790 

C      FINAL  ROW  AND  COLUMN  INTERCHANGE  011800 

K  =  N  011.810 

100  K=(K-1)  011820 

IF(K)  150,150.103  011830 

103  I=LIK)  011840 

IF(I-K)  120,120.105  011850 

105  DO110  J  =  1,N  011860 

HOLD=A(J,K)  011870 

A( J,K)=-A( J.I )  011880 

110  A( J,I)=HOLD  011890 

120  J=M(K)  011900 

IF(J-K)  100,100,125  011910 

125  D0130  1=1, N  011920 

HOLD=A(K,I)  011930 

A(K.I)=-A( J, I )  011940 

130  A{J,I)=HOLD  011950 

GO  TO  100  011960 

150  RETURN  011970 

END  011980 

C  011990 

C  012000 

SURROUTINF  HANGER  (CLX  ,CLY »CT7 .U )  0*2010 

DIMENSION  U(6.6)  012020 

DO  1  J=l,6  012030 

DO  1  K=l,6  012040 

1  U( J,K)=0.0  012050 

U(l»l)=1.0  012060 

U<2,2)=1.0  012070 

U(3,3)=1.0  012080 

U<4,4)=1.0  012090 

U(5,5)=1.0  012100 

U(6»6)=1.0  012110 

U(4.3)=CTZ  012120 

U(5.2)=-CLY  012130 

U(6»1)=CLX  012140 

RETURN  012150 

END  012160 

C  012170 

C  012180 

SUBROUTINE  HANGEO  (CTX .CLZ ,CTY ,U )  012190 

DIMENSION  U(6»6)  012200 

DO  1  J=l,6  012210 

DO  1  K=l,6  012220 

1  U(J.K)=0.0  012230 

U(ltl)«1.0  012240 

U«2,1)=CTX  012250 

U(2,2)=1.0  012260 

U(3,3)=1.0  012270 

U(4,4)=1.0  012280 

U(5,4)=CTY  012290 

U(5,5)=1.0  012300 

U(6,3)=-CLZ  012310 


m 


03 


U(6.6)=1.0 

RETURN 

END 


SUBROUTINE  CFIELO  ( A JT .R »Z ,6»SDR I t RHO.THETA, D.OI . FFAC. A ) 
DIMENSION   A(6.6) 


PHI=THETA/Z 

RHOD=RHO 

IF(RHO-0.0)11.11»12 

11  V=-1.0 
GO  TO  1 

12  V=1.0 

1  CP=COSF(PHI ) 
SP=SINF(PHI ) 
W=1.0/(  G  *AJT) 
F1=(PHI»CP+SP) /2.0 
F3=(SP-CP*PHI 1/2.0 
F5=(2.0-2.0*CP-PHI*SP 1/2.0 
F6=(2.0*PHI+PHI*CP-3.0*SP>/2.0 
RHO=ABSF(RHO) 

DO  2  J=1.6 
DO  2  K  =  1.6 

2  A(J»K)=0.0 
IF(SDRI >13tl3»14 

13  AREA=< 3.141592654/4.0 )*(D*»2-DI**2)*FFAC 
SD=RHO»THETA/(G*Z*AREA » 

GO  TO  15 
SD=0.0 


14 
15 


A(l 
A(l 
A(l 
All 
Ad 
A(2 
A(2 
A(2 
A(3 
A(3 
A(3 
A(3 
A(3 
A(3 
A(4 
A(4 
A(4 
A(4 
A(4 
A(5 
A(5 
A(5 
A(6 


RHO=RHOD 

RETURN 

END 


=  CP 

=<W*RH0»F1 )-<RHO*F3»R) 

=  V*SP 

=V*(W+R)*(RHO*PHI*SP)/2.0 

«V*(W+R)*(RHO**2)*F3 

=  CP 

*V*SP 

»V*RHO*( l.O-CP) 

»-A(2.6) 

»-1.0*A(l»6) 

■  1.0 

=RHO*SP 

-(  (R»(RHO**2  )*PHI*SP)/2.0):-W*RHO**2*F5 

*(R*RHO«*3*F3)-<W*F6*RHO**3)-SD 

=-V*SP 

=-V*(W+R)*RHO*PHI*SP/2.0 

=  CP 

=(R*RHO*Fl)-(W*RHO*F3) 

=A(3.5) 

=-V*SP 

»CP 

=A(3»4) 

=  1.0 


012320 
012330 
012340 
012350 
012360 
012370 
012  380 
012390 
012400 
012410 
012420 
012430 
012440 
012450 
012460 
012470 
012480 
012490 
012500 
012510 
012520 
012530 
012540 
012550 
012560 
012570 
012580 
012590 
012600 
012610 
012620 
012630 
012640 
012650 
012660 
012670 
012680 
012690 
012700 
012710 
012720 
012730 
012740 
012750 
012760 
012770 
012780 
012790 
012800 
012810 
012820 
012830 
012840 
012850 
012860 
012870 


04- 


C  012880 

SUBROUTINE  POINO  ( AMU  ,  AI X ,AI Y , OMEGA, TL ,Z .SDR  I ,B )  012890 

DIMENSION  BI6.6)  012900 

AM=AMU*TL/(Z-1.0)  012910 

DO  1  J=1.6  012920 

DO  1  K  =  l»6  012930 

1  B(JtK)=0.0  012940 

B(l,l)=1.0  012950 

B(2»1)=-AM»( AIX*0MEGA)*»2  012960 

B(2.2)=1.0  012970 

B<3,3)=1.0  012980 

B<4,4)=1.0  012990 

IF(SDRI-0.)9.9.10  013000 

9  B(5.4)=-AM*( AIY*0MEGA)«#2  013010 

GO  TO  11  013020 

10  B<5,4)=0  013030 

11  B(5,5)=1.0  013040 
B(6«3)=AM*0MEGA**2  013050 
B(6,6)=1.0  013060 
RETURN  013070 
END  013080 

C  013090 

C  013100 

SUBROUTINE  CFIELD  (  R »Z  ,G»SD»RHO»THETA.D.DI .E .A )  013110 

DIMENSION  A(6.6)  013120 

PHI=THETA/Z  013130 

RD=ABSE(RHO)  013140 

IF(RHO)ltlt2  013150 

1  V=-1.0  013160 
GO  TO  3  013170 

2  V=1.0         •  013180 

3  CP=COSF(PHI )  013190 
SP=SINF(PHI)  013200 
DO  4  J=1.6  013210 
DO  4  K=l,6  013220 

4  A(J.K)=0.0  013230 
A(1.1)=CP  013240 
A(1.2)=-V*SP  0L3250 
A(1.3)=-V*RD*( 1.0-CP)  013260 
A(1,4)=-V*RD**2*R*(PHI-SP)  013  270 
AR=3.141592654*(D*»2-DI**2)/4.0  013280 
W=1.0/(E*AR)  013290 
F3=(SP-PHI*CP)*.5  013300 
F1=(PHI*CP+SP)*«5  013310 
F5=(2.0-2.0*CP-PHI*SP)*.5  013320 
F6=(2.0*PHI+PHI*CP-3.0*SP)*.5  013330 
A( lt5)=V*(RD*W»PHI*SP*.5-R*F5»RD**3)  013  340 
A(  1  ,6)=RD*W*F1  +  R*F6*RD**3  013350 
A(2.1)=V*SP  013360 
A(2.2)=CP  013370 
A(2.3)*RD*SP  013380 
A(2.4)=R*RD**2*(1.0-CP)  013390 
A(2»5)=RD*F3*(W+R*RD**2)  013400 
A«2,6)=A( 1,5)  013410 
A(3,3)=1.0  013420 
A(3,4)=RD*R*PHI  013430 


OS 


A(3,5)=A(2»4) 

A(3.6)=A( 1,4) 

A<4,4)=1.0 

A(4,5)=A< 2,3) 

A(4,6)=A( 1,3) 

A(5.5)=CP 

A(5,6)=A( 1,2) 

A(6»5)=A<2,1 ) 

A(6»6)=»CP 

RETURN 

END 

c 

c 

013440 
013450 
013460 
013470 
013480 
013490 
013500 
013510 
013520 
013530 
013540 
013550 
013560 
SUBROUTINE  POINT  ( AMU  ,  A  I Y , OMEGA ,TL,Z ,R I ,B )  013570 

DIMENSION  8(6,6)  013580 

AM=AMU*TL/(Z-1. )  013590 

DO  1  J=l,6  013600 

DO  1  K=l»6  013610 

1  B(J,K)=0.0  013620 
B(l,l)=1.0  013630 
B(2,2)=1.0  013640 
B(3,3)=1.0  013650 
IF(RI)2,2,3  013660 

2  B(4,3)=-AM*(AIY*0MEGA)»*2  013670 

3  B(4,4)=l,0  013680 
B(5.2)=AM*OMEGA**2  013690 
B<5,5)=1.0  013700 
B(6,l)=-B(5,2)  013710 
B(6,6)=1.0  013720 
RETURN  013730 
END  013740 

C  '  013750 

C  013760 

C  013770 

C  013780 

C  013790 

SUBROUTINE  SFIELD  ( D I , TL,R ,Z ,G ,SD,D, E, FFAC, A )  013800 

DIMENSION  A(6,6)  013810 

ARA=( 3. 141 592654/4.0 )*(D**2-D 1**2)  013820 

AREA=ARA*FFAC  013830 

DL  =  TL/Z  013840 

DO  1  J»2,5  013850 

1  A(1.J)»0.0  013860 
A(ltl)»1.0  013870 
A(1»6)*DL/(E*ARA)  013880 
A(2,l)*0.0  013890 
A(2,2)=1.0  013900 
A(2,3)»DL  013910 
A(2»4)-DL**2*R/2.0  013920 
IF(SD)2,2,3                                                           013930 

2  A(2,5)=DL**3*R/6.0-DL/(G*AREA)  013940 
GO  TO  4  013950 

3  A(2»5)*DL**3*R/6.0  013960 

4  A(2,6)«0.0  013970 
A(3,l)=0#0  013980 
A(3t2)»0.0  013990 


06 


A(3.3)=1.0 

A(3.4)=DL»R  ■  • 

A(3.5)=A<2,4) 
A<3.6)=0.0 
DO  5  J=l»6 
A(4.J)=0.0 
A<  5»J>=0.0 
>  A(6»J)=0.0 
A<4.4)=1.0 
A(4,5)=DL 
A(5t5)«1.0 
A(6»6)=1.0 
RETURN 
END 

SUBROUTINE  SFIELO  <D  I  .  TL  ,AJT ,R ,Z.G ,SDR I .D.FFAC. A) 

DIMENSION  A(6.6)  Kijlii,llCi:,r 

AREA=( 3.141592654/4.0 >»<D**2-DI**2)*FFAC 

DL=TL/Z 
DO  2  J  =  l»6 
DO  2  K  =  l»6 
2  A<J.K)»0.0 
A(l»l)=1.0 
A(2»2)«1.0 
A(3t3)=1.0 
A<4.4)*1.0 
A(5.5)=1.0 
A(6,6)=1.0 
IF(SDRI-0.)7t7.8 

7  A(3»6)*DL**3*R/6.0-DL/(G*AREA) 

GO  TO  1 

8  A(3.6)=DL«*3*R/6.0 
1  A(1.2)=DL/(G*AJT) 

A( 3.4)=DL 
A(3.5)=DL**2»R/2.0 

A(4.5)=DL*R 

A(4»6)=A<3.5) 

A(5.6)=DL 

RETURN 

END 

SUBROUTINE  MATMUL  IA»B.U) 

DIMENSION  A(6.6>.B(6.6).C<6.6),U<6»6) 

DO  1  J  =  l»6 
DO  1  K-l»6 
C(J.  0=0.0 
DO  1  L=l»6 

1  C(  J.K)=C<  J.O+BI  J.L)*A(L»K> 
DO  2  J=1.6 

DO  2  K=l«6 

2  U( J»K)=C( J.K) 
RETURN 

END 


014000 

014010 

014020 

014030 

014040 

014  050 

014060 

014070 

014080 

014090 

014100 

014110 

014120 

014130 

014140 

014150 

014160 

014170 

014180 

014190 

014200 

014210 

014220 

014230 

014240 

014250 

014260 

014270 

014280 

014290 

014300 

014310 

014320 

014330 

014340 

014350 

014360 

014370 

014380 

014  390 

014400 

014410 

014420 

014430 

014440 

014450 

014460 

014470 

014480 

014490 

014500 

014510 

014520 

014530 

014540 

014550 


107 


c  014560 

SUBROUTINE  F  I  NBR  A(  U.VV. V.MAT  »  A  ,Z  )  014570 

DIMENSION  U(6»6) . A < 6 .6 ) » VV < 6 .3 ) »V(6.3)  014580 

TYPE  DOUBLE  VV.V  014590 

IF(MAT)  65.64.65  014600 

64  DO  67  J  =  l,6  014-610 
DO  67  K=1.3  014620 
V(J.K)=0.  014630 
DO  67  L  =  1.6  014640 

67  V( J.K)=V< J,K)+U(J.L)*VV<L.K>  014650 

RETURN  014660 

65  N=Z-1.  014670 
DO  5  M=1.N  014680 
DO  4  J=l»6  014690 
DO  4  K=1.3  014700 
V(J»K)*0.  014710 
DO  4  L=1.6  014720 

4  VI J.K)=V< J.K)+U< J.L)*VV(L.K)  014730 
DO  5  J=l«6  014740 
DO  5  K  =  1.3  014750 

5  VV( J»K)=V( J.K)  014760 
DO  7  J=1.6  014770 
DO  7  K=l»3  014780 
V(J.K)=0.  014790 
DO  7  L=1.6  014800 

7  V( J»K)=V( J.K)+A( J.L)*VV(L.K)  014810 

RETURN  014820 

END  014830 

C  014840 

C  014850 

SUBROUTINE  FINMAT  (U.UU. V.MAT , A.Z )                                     014860 

DIMENSION  U(6. 6)  .  A( 6 . 6) . V( 6.3 ) »UU ( 6 .3 )                                014870 

TYPE  DOUBLE  UU.V  014880 

IFIMAT)  65,64*65  014890 

64  DO  67  J=l,6  014900 
DO  67  K=l,3  014910 
V(J.K)*0.  014920 
DO  67  L=1.6  014930 

67  V( J.K)=V( J»K)+U( J.L)*UU(L.K)  014940 

RETURN  014950 

65  N=Z-1.  014960 
DO  5  M»1,N  014970 
DO  4  J=l,6  014980 
DO  4  <=1.3  014990 
V(J.K)«0.  015000 
DO  4  L=1.6  015010 

4  V(  J.K)=V(  J.O+U(  J.L)*UU<L.IC)  015020 
DO  5  J=1.6  015030 
DO  5  K=l#3  015040 

5  UU( J»K)*V( J.K)  015050 
DO  7  J=1.6  015060 
DO  7  K*1.3  015070 
V(J.K)=0.  015080 
DO  7  L=1.6  015090 

7  V( J.K)*V( J.K)+A( J.L)*UU(L»<)  015100 

RETURN  015110 


08 


END 


SUBROUTINE  BRANCH  ( R3 . VV»PHI ,U ) 

DIMENSION  R3<3.3).VV<6.3).U(6,6)»R1(25»25>»LL<25)«M 
1R(3.3) »G1 (3.3 ) .G2< 3.3) 
TYPE  DOUBLE  VV 
DO  4  L=1.3 
Rl( l.L)=VV( l.L) ■ 
Rl  (2»L)=VV(2.L) 
R1(3.L)=VV(3.L) 
R2( l.L)=VV(4.L) 
R2(2»L)=W(5.L) 

4  R2(3.L)=VV(6.L) 
SP=SINF(PHI ) 
CP=COSF<PHI ) 
G1(1.1)=CP 
GK1.2)  =-SP 
Gl(1.3)=0. 
G1(2.1)=SP 
G1(2.2)=CP 
Gl(2»3)=0. 
GK3.1  J=0. 
Gl(3.2)=0. 
Gl(3.3)=l. 

G2( 1 • 1 >  =1  . 

G2( 1.2)=0. 

G2(1.3)=0. 

G2(2.1>=0. 

G2(2»2)=CP 

G2(2»3)=-SP 

G2(3,l)=0. 

G2(3.2)=SP 

G2(3.3)=CP 

CALL  INVERT  ( Rl .3,0, LL .MM ) 

DO  5  J=l,3 

DO  5  K=l,3 

R3( J,K)=0.0 

DO  5  L=l,3 

5  R3(J»K)=R3(J,K)+R1(J,L)*G1(L,<) 
DO  6  J=1.3 

DO  6  K=l»3 
Rl ( J.K)=0.0 
DO  6  L=l»3 

6  Rl(J»K)=RKJ,K)+G2tJ»L)*R2(L,K) 
DO  7  J=l,3 

DO  7  K=l,3 
R( J»K)=0.0 
DO  7  L=l,3 

7  R(  J.K)=R( J.KI+R1I J,L)*R3(L.<) 
DO  8  J*l,6 

DO  8  K=l»6 

8  U(J»K)=0.0 
U(l.l)=1.0 
U(2.2)=1.0 
U(3,3)=1.0 


I  (  2  5  )  .  R  2  (  3  .  3  )  * 


015120 
015130 
015140 
015150 
015160 
015170 
015180 
015190 
0152  00 
015210 
015220 
015230 
015240 
015250 
015260 
015270 
015280 
015290 
015300 
015310 
015320 
015330 
015340 
015350 
015360 
015370 
015380 
015390 
015400 
015410 
015420 
015430 
015440 
015450 
015460 
015470 
015480 
015490 
015500 
015510 
015520 
015530 
015540 
015550 
015560 
015570 
015580 
015590 
015600 
015610 
015620 
015630 
015640 
015650 
015660 
015670 


0? 


U(4.4)=1.0 

U(5.5)=1.0 

U(6»6)=1.0 

DO  9  L=l»3 

U(4,L)=R( l.L) 

U(5.L)=R(2.L) 

9  U(6.L)=R(3tL) 

RETURN 

END 

C 

C 

SUBROUTINE  BRANCO  < R 3 . VV.PHI »U ) 

DIMENSION  R3(3*3).U(6.6).VV(6»3>»R1I25.25) 

*R2(3.3)  t 

TYPE  DOUBLE  VV 

1G1(3.3).G2(3.31.S(3.3).LL<25>,MM<25) 

DO  4  K=1.3 

Rl  (  l.K.)=VV<  1  tK) 

R1(2.K)«VV(3.K) 

R1(3.K)=VV(4,K) 

R2(1»K)=VV(2.K) 

R2(2#K)=W(5.IC) 

4  R2(3.K)=VV(6.K) 

SP=SINF(PHI ) 

CP=COSF(PHI ) 

G1(1.1)=CP 

Gl<1.2)=0. 

G1(1»3)*-SP 

Gl(2.1)=0. 

Gl(2.2)=l. 

Gl<2»3)«0. 

G1(3.1)=SP 

Gl(3.2)=0. 

G1(3.3)=CP 

G2(1.1)=CP 

G2(1»2)=SP 

G2(l»3)-0. 

G2(2.1)=-SP 

G2(2»2)=CP 

G2(2*3)>0. 

G2(3.1)=0. 

G2<3,2)=0. 

G2(3»3)=l. 

CALL  INVERT  ( Rl , 3 t D. LL .MM ) 

DO  5  J*1.3 

DO  5  K=1.3 

R3( J»K)=0.0 

DO  5  L=1.3 

5  R3<  J.K)=R3<  J.K1+RK  J.L)»G1(L.K) 

DO  6  J=1.3 

DO  .6  K  =  l»3 

R1(J.K)=0.0 

DO  6  L=1.3 

6  R1(J.K)=R1(J»K)+G2(J.L>*R2(L»K) 

DO  7  J=l»3 

DO  7  K=l»3 

S( J.K)=0.0 

015680 
015690 
015700 
015710 
015720 
015730 
015740 
015750 
015760 
015770 
015780 
015790 
015800 
015810 
015820 
015830 
015840 
015850 
015860 
015870 
015880 
015890 
015900 
015910 
015920 
015930 
015940 
015950 
015960 
015970 
015980 
015990 
016000 
016010 
016020 
016030 
016040 
016050 
016060 
016070 
016080 
016090 
016100 
016110 
016120 
016130 
016140 
016150 
016160 
016170 
016180 
016190 
016200 
016210 
016220 
016230 


!I0 


40 


43 


44 


42 


41 


15 
14 

10 


DO  7  L=lt3 
S«J.K)=S(J.<)+RHJ«L)*R3(L.K) 

DO  8  J=1.6 

DO  8  K=l»6 

U( J.K)=0.0 

U(  1.1) =1.0 

U(2.1)=S(1.1) 

U(2.2)=1.0 

U(2.3)=S<1.2) 

U(2.4)=S( 1.3) 

U(3.3)=1.0 

U(4.4)=1.0 

U(5.1)=S(2.1> 

U(5.3)=S<2.2) 

U(5.4)=S(2.3) 

U(5.5)=1.0 

U(6.1)=S(3.1  ) 

U(6.3)=S(3.2) 

U(6.4)=S<3.3> 

U(6.6)=1.0 

RETURN 

END 


SUBROUTINE  STATEMI SV I .UU.D1 .D2 .PI .P2 .PP .PM.SS ) 

DIMENSION  SVI (6.1 ) »UU(6.3) 

TYPE  DOUBLE  UU.P1 . P2 »D1 .02 

M  =  0 

DO  1  J=l»6 

DO  1  K=l»3 

UU( J.<)=0. 

IF(PM)  40.40.41 
DO  42  J  =  l»6 

IF(SVKJ.D)  42.42.43 

M  =  M+1 

UU( J»M)=1. 

IF(M-l)  42.44.42 

UU( J.2)=-D1 

UU(J.3)=-D2 

CONTINUE 

RETURN 

P1=P1+D1 

P2«P2+D2 

IF(SS)  14.15.15 

IF(SS-1.)  12.11.10 
IF(PP-1.)  12.11.10 

DO  4  J=1.6 
IF(SVKJ.D)  4.4.2 

2  M=M+1 
UU(JtM)«l« 
IFJM-2)  4.3.5 

3  UU(J.D=P1 
GO  TO  4 

5  UU(J.D=P2 

4  CONTINUE 


016240 

016250 

016260 

016270 

016280 

016290 

016300 

016310 

016320 

016330 

016340 

016350 

016360 

016370 

016380 

016390 

016400 

016410 

016420 

016430 

016440 

016450 

016460 

016470 

016480 

016490 

016500 

016510 

016520 

016530 

016540 

016550 

016560 

016570 

016580 

016590 

016600 

016610 

016620 

016630 

016640 

016650 

016660 

016670 

016680 

016690 

016700 

016710 

016720 

016730 

016740 

016750 

016760 

016770 

016780 


I  I 


11 

22 


26 
25 
?■* 
21 

1? 

32 


36 

35 
3^ 
31 


.DO  2  1  J  =  l»6 
IF( SVI ( J»l )  )  21.21  .22 
W=M+1 

UU( JtM)=l. 
IF(M-l)  2  5,26.25 
Ul.'(  J.3)  =P1 
IF(M-2J  21*23*21 

'!!'(  J,T  )  =D? 

roMT I NUF 

RFTURN 

DO  31  J =1.6 

[F(5Vi( J»l) 1  31.31  .32 

V=M+1 

UU(J»M1 =1. 
IF(M-l)  3  5.36.3  5 
UU( J.  2) =P2 
IF<m_? )  3  1,^1,33 
UU( J. 2 )=P1 
CONT INUF 
RFfu"N 


?2 

?1 
2  3 

1C 

11 
12 

3" 


31 


32 


^UR 

DIM 

TYP 

Sl  = 

B2  = 

B?  = 

01  = 

D?  = 

IF( 

1F( 

Ic( 

DO 

00 

V(  J 

IF( 

1F( 

IF( 

Bl  = 

92  = 

00 

R]  = 

R?  = 
r:,= 
GO 
Bl  = 
B2  = 
B3  = 
01  = 


ROUTIN 

ENSICN 

F  DOl'B 

0. 

0. 

0. 

0. 

0. 

PM)  3  3 

PRO!  ? 

L-l  )  ? 

1  0  *  J  =  1 

m  <=i 

»K)=YY 

SS)  11 

pp-1. ) 

SS-1.1 

V(  1  ,  1  ) 
V  (  2  .  1  ) 
vn,  1) 
TO  3 1 
Rl +V( 1 
=?+V(? 
R3  +  V< * 
TO  33 
31+V( 1 
B2+V(2 
•^3  +  v  (  3 
0 1  *  R  1 
?.1*02  + 
«1  *-Q?  + 


F  ^TATFB(Vy«»VV»N,YY,L»B3.C,Dl.n2,PP,PM,S-»Ql,Q2.Q3»FINK) 
VV(6.3).YY(3,3.25),V(3.3)  ..SV6I6.22  ) 

LF  VV 


.40,22 

1  .20,21 

3,20,2^ 

,^ 

,^ 

(  J  .  <  .  N  ) 

,12.12 

32.31,30 

33.31,30 
+  V(  1  ,2 )*D1+V( 1  ,^)*D2  +  01 
+V(2»2)*D1+V(2*3)*D2+P2 
+V<3»?J*D1+V(  3»'»)*D2+p^ 

,i  )*D?+V<  !  ,->  )«0!  +V(  1  ,  ?) 
,  1  )*r'?  +  v(  ?,  ■'  )  *D1+V(2  ,  ?  ) 
,  1  )  *D2+V(  ?,-*)  *01+V(  ?,?  ) 

. 1 I*D1+V( 1,2) *D2+V( 1,3) 
,  1  lOl+VI  2,2)  »D2+V12»3  ) 
,  1  )*Dl  +  v<  ?,  2  )  *r;?  +  V(  3 ,  i) 

c2 

n  -i 


016790 
016800 
016810 
016820 
016830 
016840 
016850 
016860 
016870 
nif-BBO 
016890 
016900 
016910 
016920 
016930 
016940 
016950 
016960 
016970 
016980 
016990 
017nno 
017010 
017020 
017030 
017040 
017050 
017060 
017070 
017080 
017090 
017100 
017110 
017120 
O17130 
017140 
017150 
017160 
017170 
017180 
017190 
017200 
017210 
017220 
017230 
017240 
017250 
017260 
017270 
017280 
017290 
017300 
017310 
017320 
017330 


I  12 


20 

1 


2 

3 

35 

34 
36 

7 
8 


9 

4 

40 
54 
52 

55 

51 

41 

43 

44 
42 


20 
10 

12 

11 


M=l 

DO  1  J=l»6 

DO  1  K= 1.3 

VV< J.K)=0.0 

DO  4  J=l,6 

IF(SVB( J.N) )  4,4.2 

IF(M-1)7,3,7 

IF(PM)  34.35.35 

VV( J,1)*Q1 

GO  TO  36 

VV< J.l )=1. 

M=M+1 

GO  TO  4 

IF(M-2)  4,8,9 

VVIJ.l )=Q2 

VV(J»M)«1. 

M=M+1 

GO  TO  4 

VV( J.l )=Q3 

VV(J.M)=1. 

CONTINUE 

RETURN 

IF(FINK)  54,51.54 

IF(BRC)  51.51.52 

DO  53  J=1.3 

DO  53  K»l»3 

V<  J.K)«YY( J.K.N) 

D1=(V<1.2)/V( 1.1)+V(2.2)/V(2»1)+V(3.2)/V<3.1 ) )/3. 

D2=(V(1.3)/V(1.1)+V(2.3)/V(2.1)+V(3.3)/V(3.1 ) )/3. 

M=0 

DO  41  J=l,6 

DO  41  K  =  1.3 

VV< J.K)=0. 

DO  42  J*1.6 

IF(SVB( J.N) )  42.42.43 

M*M+1 

VV( J.M)=1, 

IF(M-l)  42.44.42 

VV( J.2)=-D1 

VV( J.3)=-D2 

CONTINUE 

RETURN 

END 


SUBROUTINE  BRCORJ R3.UU.XX. JK.VV.PM ) 

DIMENSION  R3(3.3),UU<6,3)»XX(3.3) .UUU«3.3) »VV(6,3) 

TYPE  DOUBLE  VV.UU 

IF(PM)  20.21.20 

IF(JK)  10.10.11 

DO  12  J=1.3 

DO  12  K=l,3 

UUU( J.K)=UU( J.K) 

GO  TO  3 

DO  13  K=1.3 

UUU( l.K)=UU< l.K) 


017.340 
017350 
017360 
017370 
017380 
017390 
017400 
017410 
017420 
017430 
017440 
017450 
017460 
017470 
017480 
017490 
017500 
017510 
017520 
017530 
017540 
017550 
017560 
017570 
017580 
017590 
017600 
017610 
017620 
017630 
017640 
017650 
017660 
017670 
017680 
017690 
017700 
017710 
017720 
017730 
017740 
017750 
017760 
017770 
017780 
017790 
017800 
017810 
017820 
017830 
017840 
017850 
017860 
017870 
017880 
017890 


I  13 


UUU(2.K)=UU<3»K)  017900 

13  UUU(3.K)=UU(4,K)  017910 

3  DO  1  J=1.3  017920 

DO  1  K=l,3  017930 

XX<J.K)«0.  017940 
DO  1  L*l»3                                                            .  017950 

1  XX( J,K>-XX< J»K)+R3< J.L)*UUU(L,K)  017960 
RETURN  017970 

21  IF'(JK)  25,25.22  017980 
25  DO  23  J  =  1.3  017990 

DO  23  K»l,3  018000 

23  XX( J,K)«VV( J»K)  018010 
RETURN  018020 

22  DO  24  K«l»3  018030 
XX(1»K)«VV<  1>»K)  018040 
XX«2,IC)»VV{3»tC)  018050 

24  XX(3.K)=VV(4,K)  018060 
RETURN  018070 
END  018080 

C  018090 

C  018100 
SUBROUTINE  DELMA( UU, SVE tV.X2 » Y2 ,KKK,X1 ,FF ,PP .HH.RR ,SS . Yl ,REM,F I NK.   018  110 

1PM.D1.D2)  018120 

DIMENSION  UU<6»3) »V<6»3) »SVE<6»1> »D(3»3)  018130 

TYPE  DOUBLE  UU» V.D.X 1 »X2 » Yl . Y2 .DV.D1 »D2  018140 

L=0  018150 

X1*0.  018160 

X2-0.  018170 

Y2  =  0.  018180 

Y1*0.  018190 

Dl-0.  018200 

D2=0.  018210 

DO  111  J*l»3  018220 

DO  111  K«l»3  018230 

1U  D(J,K)«0.  018240 

DO  91  N-1.6  018250 

IFISVE(N.l))  91,2.91  018260 

2  L«L+1  018270 
DO  92  K*l,3  018280 

92  V(L.K)*UU(N.K)  018290 

91  CONTINUE  018300 

IF(PM)  66,64*61  018310 

64  IF(FINK)  63,62,63  018320 

62  D1=(V( 1,2)/V(1,1)+V(2,2)/V(2,1)+V(3,2)/V(3,1) )/3.  018330 
D2«(V(l,3)/V(l,l)+V(2,3)/V(2,l)+V(3,3)/V(3»l))/3.  018  340 
FINK-1.  018350 
RETURN  018360 

63  FINK=0.  018370 
RETURN  018380 

66  RETURN  018390 

61  IF(SS)  311,312,312  018400 

311  IF(PP-1.)400, 401,402  018410 

312  IFJSS-1.)  313,314,315  018420 
402  IF(RR-1.)  300,403,300  018430 
315  IF(RR-1.)  340,341,340  018440 
300  IF(V(3,3))  93,94,93  018450 


340  IF(V<3.2))  93.94,93  018460 

93  DO  95  J=2t3  018470 
DO  95  K=»1.3  018480 

95  D( J.K)=V< J»K)  018490 
DV=V(1.2)*D(2»3)-V(1 t3)*D(2.2)  0185  00 
IF<DV)  98,94.98  018510 

94  RR=1.0  018520 
HH=1.  018530 
RETURN  018540 

403  IF(V(2.3))  96.200.96  018550 

341  IF(V(2.2) )96.200.96  018560 

96  DO  97  K=1.3  018570 
D(2.K)=V(3.K>  018580 

97  D(3.K)»V(2fK,)  018590 
DV=V(1.2)*D(2.3)-V(1.3)*D<2.2)  018600 
IF(DV)  98.200.98  018610 

98  X1=X1-(V( 1.1)*D(2.3)-V(1.3)*D(2.1) )/DV  018620 
X2=X2-(V(1.2)*D(2.1)-V{1.1)*D(2.2) )/DV  018630 
IF(X1)  100.102.100  018640 

102  IFCX2)  100.101.100  018650 

101  FF*1.0  018660 

RETURN  018670 

100  IF(SS)  354,355.355  018680 

354  Y2=Y2-D(3.1)/D(3.3)-X1*D(3.2)/D(3,3)  018690 
RETURN  018700 

355  Y1=Y1-D<3.1 )/D(3.2)-X2*D(3.3)/D(3.2)  018710 
.  RETURN  018720 

99  IF(SS)  20,21.21  018730 

20  PP=1.  018740 
GO  TO  23  018750 

21  SS=1.  018760 
23  RR=0.  018770 

HH«1.  018780 

RETURN  018790 

314  IF(RR-1.)  330.331.330  018800 

401  IF(RR-l.O)  301.404.301  018810 

330  IF(V<3,1)>  81.82,81  018820 
301  IF(V<3.2))  81.82.81  018830 

81  DO  83  J=2.3  018840 
DO  83  K  =  l,3  018850 

83  D(J,K)=V< J.K)  018860 
DV=V(1, 1)*D(2»2)-V(1,2)*D(2,1)  018870 
IF(DV)  86,82,86  018880 

82  RR=1.0  018890 
HH«1.  018900 

RETURN  018910 

331  IF(V(2.1))  84,70,84  018920 

404  IF(V(2,21>  84,70,84  018930 

84  DO  85  K=l,3  018940 
D(2.K)»V(3,K)  018950 

85  D(3.K)=V(2.K)  018960 
DV=V(1.1)*D(2.2)-V(1.2)*D(2.1)  018970 
IF(DV)  86.70.86  018980 

86  X1=X1-(V( 1.3)*D(2.2)-V(1.2)*D(2»3) )/DV  018990 
X2=X2-(V<1.1)*D<2.3)-V(1.3)*D(2»1) )/DV  019000 
IF(X1)2C3,201.203  019010 
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201  IF(X2)  203t2O2,2O3  019020 

202  FF=1.0  019030 
RETURN  019040 

203  IF(SS) 352,353,353  019050 
353  Y1=Y1-D(3. 31/0(3.1 )-X2*D( 3,2)/D(3.1 >  019060 

RETURN  019070 

352  Y2=Y2-D(3»3)/D(3»2 >-Xl«D( 3,1 )/D< 3,2)  019Q80 

RETURN  019090 

70  IF(SS)  25.24.24  019100 

24  SS=2.  019110 
GO  TO  26  019120 

25  PP=2.  019130 

26  RR  =  0.  019140 
HH=1.  019150 
RETURN  019160 

313  IFIRR-1.)  320.321.320  019170 

400  IFIRR-1.)  302.405,302  019180 

320  IFIVI3.3))  71.72.71  019190 
302  IFIVI3.1) 171,72.71  019200 

71  DO  73  J=2.3  019210 
DO  73  K=1.3  019220 

73  DIJ»K)=V(J.K)  019230 
DV=VI1.1)*D(2.3)-V( 1.3)»DI2.1 )  019240 
IFIDV)  76.72.76  019250 

72  RR=1.0  019260 
HH=1.  019270 
RETURN  019280 

405  IFIVI2.11)  74,99,74  019290 

321  IFIVI2,3>)  74,99,74  019300 
200  KKK»1  019310 

RR=0.  019320 

RETURN  019330 

74  DO  75  K=l,3  019340 
D(2,K)*V(3,K)  019350 

75  D<3,K)-V(2,IC)  019360 
DV  =  VI 1 ,1 )*D(2»3 >-V(  1  ,3)*D(2,1  )  019370 
IFIDV)  76,99,76  019380 

76  X2=X2-IV<1,2)*D(2,3)-VI1.3)*D<2»2) >/DV  019390 
X1=X1-(VI1,1)*DI2,2)-D(2,1)»VI1,2))/DV  019400 
IFIX1)  213,211,213  019410 

211  IFIX2)  213,212,213  019420 

212  FF=1.0  019430 
RETURN  0194<»0 

213  IFISS)  350,351,351  019450 
351  Y1=Y1-IDI3,2)+X2*DI3,1 ) )/D(3,3)  019460 

RETURN  019470 

350  Y2=Y2-(D(3,2)+X1*D(3,3 ) 1/DI3.1)  019480 

RETURN  019490 

END  019500 
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APPENDIX  E 
ACCURACY  OF  METHODS  &  SAMPLE  PROBLEMS 
E-l  General  remarks. 

Appendix  E  contains  the  results  of  vibration  analyses  of  several 
piping  systems.   The  first  group  was  performed  to  establish  the  accuracy 
of  the  methods  and  compare  them  with  the  VIPIPE  output.  The  second  group 
was  performed  to  establish  assurance  of  solution  by  comparing  the  natural 
frequencies  got  from  various  methods. 

E-2  Method  accuracy  analyses. 
E-2-1  Straight  sections. 
End  conditions. 

System  1  -  fixed-fixed. 
System  2  -  fixed-free. 
System  3  -  fixed-propped. 
System  parameters  (for  system  1,2,3). 
length  -  200  inches, 
diameter  -  2.0  inches, 
wall  thickness  -  0.125  inches. 
Intensive  properties  (for  system  1,2,3) 
density  -  490  lbs/cu-ft 

shear  modulus  -  12x10  psi 

6 
elastic  modulus  -  30x10  psi 

Mathematical  model:   Distributed  mass  with  the  effect  of  shear 

deflection  and  rotational  inertia  neglected  to  permit  comparison  with 

classical  solution. 
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Comparison  of  results. 

As  explained  in  Appendix  B-9-4,  when  the  iteration  reaches  the 
solution  acceptability  criterion,  the  corresponding  natural  frequency 
is  determined  using  linear  interpolation.   Theoretically  exact  compari- 
son values  are  obtained  by  use  of  the  same  iteration  method  with  the 
solution  acceptability  criterion  as  0.0007  rad/sec.   The  values  obtained 
using  VIPIPE  have  a  solution  acceptability  criterion  0.0185  rad/sec. 
Comparison  is  made  by  percent  difference. 

The  M-D  method  fails  the  systems  1  through  3,  because  the  purging 
factors  turn  out  to  be  all  zeros.   So  for  the  purpose  of  testing  the 
M-D  method,  modified  System  1  is  provided  below. 

The  letters  D  or  S  in  parentheses  denote  double  precision  and  single 
precision  respectively. 

Modified  System  1: 

Add  a  very  small  curved  section  to  the  System  1.   Second  section  is 
a  mathematical  model  which  has  the  same  property  as  System  1  except  the 
system  parameters. 

Second  section  parameters: 

/included  angle  of  arc  -  45  degrees 
radius  of  curvature  -  0.01  inches 
(This  value  of  radius  of  curvature  is 
not  possible  in  an  actual  system,  since  the  radius  of  curvature  less  than 
the  radius  of  the  pipe.   In  effect  subroutine  STIFCO  or  STIFOO  will  pro- 
vide the  transfer  matrix  of  a  massless  rigid  corner). 
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Results  and  comparison.* 

a.   System  1  and  modified  System  1.**  (in-plane). 


Mode 

Freq. (rad/sec)  and  difference  in  7.. 

No. 

Comparison  values 

Delta  m.  (D) 

P-M.-A(D) 

P-Mm.-F(D) 

1 

75.09994540 

75.09994544 
0     % 

75.09994550 
0     % 

75.09994571 
0     % 

2 

207.01589137 

207.01589036 
5.0xl0~77. 

207.01589094 
0     % 

207.01589228 
4.4xl0'7  % 

3 

: 
i 
i 

405.83391916 

405.83393183 
3. 0x10 'So 

405.83403022 
2.4xl0"5% 

405.83394702 
7.0xl0"77. 

4 

j 

670.86408383 

670.86310856        670.86341454 
2.8xl0"4%          1.0xl0~57. 

670.86345124 
9.1xl0"5% 

j 

■ 
5 

1002.15520225 

1002.19978814       1002.18735784 
4;7xl0"3%          3.2xl0"3% 

1002.19300449 
3.8xlO"3% 

f 

j 

6 

1399.70436269                       1399.59981662 

7.5xlO-3% 

1399.64448571 
4.3xl0"\ 

I 
j 

■ 

1       \   1863.51172599 

J 

1959.59556134 
5.1    7. 

Mode 

Freq. (rad. sec)  and  difference  in  % 

No. 

Delta  m.  (S)      P-M  m.-A(S) 

M-D  m.(D)*** 

M-D  m.(S)*** 

1. 

75.09994535       75.09994865 
0    7.        4.0xl0"6  7. 

75.09994539 
0     7. 

75.09994534 
0      7, 
207.01588891 
0     7, 

2 

207.01589228      207.01589071 
5.0xl0"87.         0     7o 

207.01588836 
1.5xl0"6  7o 

3     405.83398145      405.83403455 
1.5xl0"67o        2.1xl0"5  7. 

405.83395151 
9.5xl0"6  7o 

405.83401120 
2.0xl0-5  7. 

4     670.86345124 
7.5xl0'\ 

670.86256288 
2.2xl0"4  7o 

670.86374475 
5.1xl0"5  7. 

670.86118054 
4. 3x10 "4  7. 

5 

1002.23499946 
8.0xl0"3  7. 

1002.19909306 
4.4xl0"4  % 

1002.21779856 
6.3xl0"3  % 

! 6 

i 

1399.17084415 
4. 5x10" 2  7o 

1399.94044241 
3. 4x10" 2  % 

♦Sources  of  comparison  frequencies  for  system  1  through  3  are  given  in  £-2-2. 
(Additional  footnotes  on  next  page). 
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**In  all  the  tables  which  follow,  simple  and  evident  abbreviations  are  used, 
For  example,  P-M  m.-A(D)  means  the  P-M  method,  sub-procedure  A  in  double 
precision  arithmetic. 

***Modified  System  1. 
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b.      System  2. (in-plane) 


Mode 
NO 


Comparison  values 


11. 80213586 


73.96272297 


207.09776595 


405.82896584 


670.86435913 


1002.1551877 


1399.70436329 


Freq.(rad/sec)    and  difference   in  % 


Delta  m.    (d) 


11.80213389 
5.5xl(T*  % 

73.96272167 
1.7xl0"6  % 


207.09776685 
4.1xl0~8  # 


405.82894132 
6.1xlO~6  # 


670.86200817 
2.2x10-5  4 


1002.26833752 
1.1x10-2  % 


1399.08849534 


2.0x10 


-2 


P-M  m.-A  (D) 


11.80214614 

9. 5x10" 5  # 

73.96272478 
1.7x10-6  # 


207.09776388 
l.lxlO"6  0 


405.82893748 
7.0xl0"6  # 


670.86326697 
1.7xlO"5  % 


1002.18809484 
4.4xlO-5  $ 


1399.61231065 
6.6xlO-5  % 


P-M  m.-P  (D) 

11.80213746 
1.4xl0-5  % 


73.96272367 
0  # 


207.09776562 
0  # 


405.82900191 
9.1xl0"6  % 


670.86325735 
1.7x10-5  £ 


1002.19715557 
4.2xl0"5  % 


1399.14599562 
4.2x10-5  # 


8 


1863.51172626 


1862.55983728 
5.1xlO"2  % 


c.   System  3. (in-plane) 


Mode 
NO 


Comparison  values 


51.75397285 


167.71502090 


349.92609071 


598.39432089 


913.12074582 


1294.10536551 


1741.34817946 


Freq.(rad/sec)    and  difference  in  # 


Delta  m.    (D) 


51.75397263 
0  # 


167.71602009 
0     % 


349.92608231 
l.lxlO"7  % 


598.39381063 
8.5x10-5  % 


913.13150640 
1.2xl0"5  % 


1293.77838653 
2.7xl0"2  # 


P-M  m.-A  (D) 


51.75397569 
0  # 


167.71602200 
0  S 


349.92610259 
3.4xKT6  % 


598.39446341 
2.5x10-5  $> 


913.13180993 
_       1.2xl0"3  $ 
1294.11223844 
5.8xl0"4  # 


1748.59011823 
4.2X10"1  % 


P-M  m.-P  (D) 


51.75397564 
0  # 


167.71602189 
0  # 


349.92610209 
3.4xlO-6  % 


598.39407344 
4.2x10-5  £ 


913.13109335 
1.2x10-5  % 


1294.11761177 
1.9xlO~4  # 
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E-2-2  Sources  of  comparison  frequencies  for  system  1  through  3. 

The  governing  fourth  order  equation  of  a  homogeneous  single  component 
straight  system 

E-X   J 
was  solved  to  give 

y  ■  a  cos  dx  +  b  coshcXx  +  c  sin  Ay-  +  &   sinh^x 
Upon  applying  boundary  conditions  for  System  1  (y  ■  y'  =  0  at  x  ■  0, 
x  -  L)  and  simplifying,  it  was  found  that  the  eigenvalues  of  the  system 
must  satisfy 

cos  of  ■  cosho<  -  1  (E-2-2-1) 


where  4/x  u)4  ^ 

^~     J    EX  (E-2-2-2) 

Upon  applying  the  boundary  conditions  for  System  2(y  ■  y*  ■  0  at  x  ■  0, 

y"  =  y"1  =  o  at  x  *  L),  it  was  found  that  the  eigenvalues  of  the  system 

must  satisfy 

cos  dt  •  cosho(  -  -1  (E-2-2-3) 

Upon  applying  the  boundary  conditions  of  System  3(y  ■  y1  ■  0  at  x  ■  0, 
y  =  y"  =  o  at  x  =  L)j  it  was  found  that  the  eigenvalues  must  satisfy 

sino(  coshO(-  cos  tf  cosh  o(   »  0         (E-2-2-4) 
The  roots  of  equation  E-2»2-l,  E-2-2-3  and  E-2-2-4,  from  which  the  natural 
frequencies  of  the  corresponding  systems  may  easily  be  found,  were  ob- 
tained through  iteration  using  a  digital  computer. 


122 


E-2-3  Curved  section  (System  4) 

System  parameters: 


radius  of  curvature 

included  angle  of  arc 

diameter 

wall  thickness 

Intensive  properties: 

density 
shear  modulus 
elastic  modulus 


50  inches 
270  degrees 
2.0  inches 
.125  inches 


556.  lbs/cu-ft 
12x106  psi 
30xl06  psi 


End  conditions: 
Mathematical  model: 

Results. 

a.  In-plane  vibration* 


fixed- fixed. 

lumped  mass  with  shear  deflection  and  rota- 
tional inertia  neglected. 


-, 

Mode 

. — — — _ — _ _ -i, 

Frequency  (rad/sec) 

No. 

Delta  m.(D) 

M-D  m.  (D) 

P-M  m.-A  (D) 

P-M  m.-C  (D) 

i-4 

66.69934902 

66.69934902 

66.69934980 

66.69934948 

2 

171.22352538 

171.22352538 

171.22352610 

171.22352592 

1   3 

335.35007083 

335.35007083 

335.35007083 

335.35007174 

4 

535.64529414  j  535.6429665 

535.64529637 

5 

773.44358487  j  773.44358842 

773.44358812 

b.  Out-of -plane  vibration* 


Mode 
No. 


i 


Frequency  (rad/sec) 


Delta  m. (D) 


35.36309735 


94.42588921 


207.89287257 


369.11199342 


M-D  m.  (D) 


35.36309735 


94.42588921 


P-M  m.  -A  (D) 


35.36309778 


94.42588653 


207.89278647 


207.89287345 


369.11278647  j  369.11278689 


571.22235666  j  571.22235577 


*See  footnote  next  page, 
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*The  fundamental  frequency  value  of  in-plane  and  out-of-plane  case  are.  . 
close  to  the  values  obtained  by  employing  equation   0=  "F(oOf(  E-VSwR 
after  obtaining  F(  o<  )  from  the  curves  in  reference   [5].   Natural 
frequency  values  from  the  three  methods  are  the  same  up  to  six  digits. 
One  can  get  more  higher  mode  frequencies  using  the  M-D  and  P-M  method. 
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E-3  Other  Sample  Problems 
A.  System  5 

System  parameters 


! 

component 

Parameter 

1 

2 

3 

4      5 

6 

7 

diameter  (in) 

2.0 

2.0 

2.0 

H     2.0 
A 

2.0 

2.0 

wall  thickness 

).75 

0.75 

0.75 

N      0.75 

0.75 

0.75 

(in) 

G 

length  (in) 

130. 

1.40 

6.0 

E     60.21 
R 

1.50 

6.28 

radius  of  curve 

4.0 

30. 

4.0 

|  (in) 

:  incl.  angle  of  arc 

.1  ... 

20. 

115. 

90. 

Intensive  properties  of  all  components  except  hangers: 

density:   556  lbs/cu-ft.* 

6 

shear  modulus:  12x10  psi 

6 
elastic  modulus:  30x10  psi 

Hanger  spring  rates: 

CLX  ,  CLY  ,  CLZ  10001b s/ in. 

CTX  ,  CTY  ,  CTZ  1000  in-lb/rad. 
End  conditions: 

both  ends  of  system  fixed. 
Mathematical  model:  Distributed  mass  with  the  effects  of  shear  deflec- 
tion and  rotational  inertia  considered. 


^ 


H.4  J' 


*This  fictitious  value  was  chosen  to  permit  comparing  with  results  of 
Fink  [4]  who  used  the  same  value. 
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Results,  (in-plane  vibration) 


Mode 
NO 

Frequency ( rad/sec) 

Delta  m.    (D) 

M-D  m.    (D) 

P-M  m.-A  (D) 

P-M  m.-C   (D) 

1 

116.89479180 

2 

304.88354398 

3 

566.43670976 

4 

826.88282292 

826.88281950 

826.88281929 

i 

5 

1126.83773518 

1126.83799669              1126.83797508 

6 

1666.58935151 

1666.58812881              1666.58878192 

^. 

7 

2099.54173011                2099.5352509 

2099.52844030 

Results  in  single  precision  arithmetic  and  comparison  with  double  precision  results, 


Mode 
NO 

Frequency^ rad/sec)    and  diff.    in  # 

Delta  m.    (s)                   M-D  m.    (b)                       P-M  m.-A  (s) 

P-M  m.-C    (S) 

1 

116.89479181 
•    0              % 

2 

304.88354658 
l.OxlO"6  % 

• 

• 

. 

3 

566.43620786 
8.9xlO-5  % 

* 

4 

826.88281494 
l.OxKT6  # 

826.88281417 
6.0x10-7  <f> 

826.88281764 
1.9x10-7    . 

5 
6 

1126.83768564 
4.4xlO=6  # 

1126.83786333 
1.2xlO-5  $> 

£__ 

1126.83745223 

4.6x10-5  # 

1666.58927426 
4.7xKT6  # 

«| 

7 

2099.53091794 
5.0xlO"3  # 

——.-.. 
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B.   System  6  (single  branch  system) 

System  parameters: 


Bra?Jch 


Intensive  properties: 


parameter 

components 

1 

2       3 

4 

dia.  (inch) 

2.0 

2.0      2.0 

2.0 

1 

wall  thickness 

(inch) 
length  (inch) 

.125 
100. 

.125     .125 
100.     .01 

.125 
100. 

incl.  angle  of 
arc.  (degree) 

0 

0       45 

0 

density 


556  lbs/cu-ft 


shear  modulus   12x10  psi 


elastic 
modulus 


30x10  psi 


End  conditions:  All  ends  fixed. 

Mathematical  model:  Distributed  mass  with  the  effects  of  shear  deflec- 
tion and  rotational  inertia  neglected. 


Result 

9 

!  Mode 

Delta  m.(D) 

Freqi 
M-D  m.  (D) 

lency  (rad/sec) 
P-M  m. -A  (D) 

No. 

P-M  m.-C  (D) 

1 

1 

194.25470558 

2 

i 

278.61078241 

3 

628.59556355 

4 

747.68138190 

747.68087208 

747.68137310 

747.68137759 

!   5 

1300.29478437 

1300.31646085 

1300.31641281 

6 

1381.79339427 

1381.86234251 

1381.86230734 

7 

1991.80356348 

1991.38946110 

1991.38946101 

8 

. . 

2257.92926800 
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